{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The autoreload extension is already loaded. To reload it, use:\n",
      "  %reload_ext autoreload\n"
     ]
    }
   ],
   "source": [
    "import torchdyn\n",
    "import torch\n",
    "import torch.nn as nn\n",
    "import matplotlib.pyplot as plt\n",
    "from torchdyn.numerics import Euler, RungeKutta4, Tsitouras45, DormandPrince45\n",
    "from torchdyn.numerics import odeint\n",
    "from torchdyn.core import ODEProblem\n",
    "\n",
    "import torchdiffeq\n",
    "import time \n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "class VanDerPol(nn.Module):\n",
    "    def __init__(self, alpha=10):\n",
    "        super().__init__()\n",
    "        self.alpha = alpha\n",
    "        self.nfe = 0\n",
    "\n",
    "    def forward(self, t, x):\n",
    "        self.nfe += 1\n",
    "        x1, x2 = x[...,:1], x[...,1:2]\n",
    "        return torch.cat([x2, self.alpha * (1 - x1**2) * x2 - x1], -1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = VanDerPol()\n",
    "x = torch.randn(1024, 2)\n",
    "t_span = torch.linspace(0, 3, 300)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Fixed--step bench"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0637514591217041\n",
      "0.06930422782897949\n"
     ]
    }
   ],
   "source": [
    "t0 = time.time()\n",
    "t_eval, sol1 = odeint(f, x, t_span, solver='rk4')\n",
    "t_end1 = time.time() - t0\n",
    "print(t_end1)\n",
    "\n",
    "t0 = time.time()\n",
    "sol2 = torchdiffeq.odeint(f, x, t_span, method='rk4')\n",
    "t_end2 = time.time() - t0\n",
    "print(t_end2)\n",
    "\n",
    "true_sol = torchdiffeq.odeint(f, x, t_span, method='dopri5', atol=1e-9, rtol=1e-9)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f7829e18400>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfEAAADbCAYAAAB9cqiZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAABWsUlEQVR4nO3dd3gUVdsG8PtsSYcQem8CUkOoAgJSRVBBpSh2ECniSxHEThMVsHyoCIqKFJWugvTeCZ2EBAIk1NASSEivu/f3x0xiEjYhwCabkOd3XXtld+bMzLOTTZ49Z86co0hCCCGEEIWPwdEBCCGEEOLeSBIXQgghCilJ4kIIIUQhJUlcCCGEKKQkiQshhBCFlCRxIYQQopAqlElcKTVXKRWmlAqw0/4sSqlj+mOVPfYphBBC5DVVGO8TV0q1BxALYAHJhnbYXyxJj/uPTAghhMg/hbImTnIngIiMy5RSDyml1iulDiuldiml6jooPCGEECJfFMokno05AP5HshmAsQBm3cW2LkqpQ0opX6XUM3kSnRBCCGFnJkcHYA9KKQ8AbQAsU0qlLXbW1z0HYLKNzS6T7KY/r0byslKqJoCtSqnjJEPyOm4hhBDifjwQSRxai8Itkj5ZV5D8C8BfOW1M8rL+86xSajuAJgAkiQshhCjQHojmdJLRAM4ppfoCgNI0zs22SikvpVRarb00gEcBnMizYIUQQgg7KZRJXCm1CMA+AA8rpUKVUm8AeAnAG0opPwCBAHrlcnf1ABzSt9sGYCpJSeJCCCEKvEJ5i5kQQgghCmlNXAghhBCSxIUQQohCq9D1Ti9dujSrV6/u6DCEEEKIfHP48OEbJMtkXV7oknj16tVx6NAhR4chhBBC5Bul1AVby6U5XQghhCikJIkLIYQQhZQkcSGEEKKQKnTXxIUQQuS/lJQUhIaGIjEx0dGhPNBcXFxQuXJlmM3mXJWXJA4Av/8O1KsHNGvm6EiEEKJACg0NRbFixVC9enVkmGhK2BFJ3Lx5E6GhoahRo0autpEkDgCvvKL9lNHrhBDCpsTEREngeUwphVKlSiE8PDzX28g1cSGEELkiCTzv3e05liQuhBCiwLt16xZmzZpll31Vr14dN27csFs5R5IkLoQQosC72ySempqah9EUHJLEM7JaHR2BEEIIG95//32EhITAx8cH7777Lt599100bNgQjRo1wpIlSwAA27dvR7t27dCzZ0/Ur18fFosFY8eORcOGDeHt7Y3vv/8+fX/ff/89mjZtikaNGiEoKAgAcPPmTTz++ONo0KABBg0ahLRZPsePH48ZM2akb/vRRx/h22+/xfbt29GhQwf06dMHdevWxUsvvYT8nhlUOrZlFBEBlC7t6CiEEKJAGzVqFI4dO2bXffr4+GRKlFlNnToVAQEBOHbsGFasWIEff/wRfn5+uHHjBlq0aIH27dsDAI4cOYKAgADUqFEDs2fPxvnz53Hs2DGYTCZERESk76906dI4cuQIZs2aha+++gq//PILJk2ahLZt22L8+PFYs2YNfv31VwDAwIED8dxzz2HUqFGwWq1YvHgxDhw4gOPHj+Po0aMIDAxExYoV8eijj2LPnj1o27atXc9NTqQmnkHypUuODkEIIcQd7N69G/3794fRaES5cuXw2GOP4eDBgwCAli1bpt+etXnzZgwZMgQmk1ZfLVmyZPo+nnvuOQBAs2bNcP78eQDAzp078fLLLwMAnnzySXh5eQHQro2XKlUKR48excaNG9GkSROUKlUq/XiVK1eGwWCAj49P+r7yi9TEM4gMCkK5Jk0cHYYQQhRoOdWYHc3d3T1X5ZydnQEARqMxV9fPBw0ahHnz5uHatWsYOHDgbfu5m33Zk9TEM4gJCXF0CEIIIWwoVqwYYmJiAADt2rXDkiVLYLFYEB4ejp07d6Jly5a3bdO1a1f89NNP6Yk1Y3O6Le3bt8eff/4JAFi3bh0iIyPT1z377LNYv349Dh48iG7dutnrbd03SeIZJFywOdObEEIIBytVqhQeffRRNGzYEPv27YO3tzcaN26MTp06Yfr06Shfvvxt2wwaNAhVq1ZNL5uWoLMzYcIE7Ny5Ew0aNMBff/2FqlWrpq9zcnJCx44d0a9fPxiNRru/v3ul8rsn3f1q3rw57TqfOAkYtO8yR7t3R5O1a+23byGEeECcPHkS9erVc3QYDmO1WtG0aVMsW7YMtWvXztNj2TrXSqnDJJtnLSs18YxfYq5fd1wcQgghCqQTJ06gVq1a6Ny5c54n8LtV5Du2MTUVaYPcmW7edGgsQgghCp769evj7Nmzjg7DpiJfE7dm6EnoGhXlwEiEEEKIu1Pkk7glJSX9uUdCggMjEUIIIe5OkU/i1gxJvHhysgMjEUIIIe6OJHG9OT0GgBspc4oLIYQoNIp8Ek9rTo/TX1Oa1IUQosBx9FSkbdq0SV/+7rvvokGDBnj33XcRHh6ORx55BE2aNMGuXbvsEt/dkN7pek083mgELBYk37oFZzc3B0clhBAio7Qk/tZbb+WqfGpqavqY6fawd+/e9Odz5sxBREQEjEYjFi9ejEaNGuGXX36x27HuRpGviac1pyfoI/Ak3GFYPiGEEPnPkVORAoCHhwcAoGfPnoiNjUWzZs0wbdo0jBs3DitXroSPjw8SEhKwceNGtG7dGk2bNkXfvn0RGxsLAFi/fj3q1q2Lpk2bYsSIEXjqqafscl7yrCaulKoCYAGAcgAIYA7Jb7OUUQC+BdADQDyA10keyauYbElrTk8ymYDkZCRlGCtXCCGEDaNGAXaeihQ+PkABnYo0o1WrVsHDwyN9KtZy5crh0KFDmDlzJm7cuIEpU6Zg8+bNcHd3x7Rp0/DNN99g3LhxePPNN7F161bUqlULzz//vN1OW17WxFMBjCFZH0ArAMOVUvWzlOkOoLb+GAxgdh7GY1Nac3qSkxMAIFFq4kIIUaDl91SkueXr64sTJ07g0UcfhY+PD+bPn48LFy4gKCgINWrUQO3ataGUSj+GPeRZTZzkVQBX9ecxSqmTACoBOJGhWC8AC6i1WfgqpUoopSro2+aLtJp4ij6dXLIM+CKEEDkrglOR5gZJdO3aFYsWLcq0/Ji9Wy0yyJdr4kqp6gCaANifZVUlAJcyvA7Vl+UbWiwAgFRJ4kIIUWA5eirS3GjVqhX27NmD4OBgAEBcXBxOnz6NunXr4vz58wjRp7vOmuTvR54ncaWUB4AVAEaRjL7HfQxWSh1SSh0KDw+3a3xpg71YXF0BAKmSxIUQosBx9FSkuVGmTBnMmzcP/fv3h7e3N1q3bo2goCC4uLhgzpw5ePLJJ9G0aVOULVv2rvabkzydilQpZQawGsAGkt/YWP8TgO0kF+mvTwHokFNzur2nIg3dsQOVO3TAhiZN0O3oURweNw7Npk2z2/6FEOJBUNSnIrWn7du346uvvsLq1attri8QU5HqPc9/BXDSVgLXrQLwqtK0AhCVn9fDgQwToOi3D1j05hohhBCioMvLwV4eBfAKgONKqWP6sg8BVAUAkj8CWAvt9rJgaLeYDcjDeGxK652OYsUAABb9nj4hhBAiL3To0AEdOnSwy77ysnf6biB9qu7syhDA8LyKITfSauKqeHEAACWJCyGEKCSK/Ihtab3TTWlJPD7ekeEIIUSBlZd9qITmbs9xkU/iab3TzW5uSAIASeJCCHEbFxcX3Lx5UxJ5HiKJmzdvwsXFJdfbFPkJUNKa041OTogHJIkLIYQNlStXRmhoKOx9m6/IzMXFBZUrV851+SKfxNM6thlMJiQqBZWY6OCIhBCi4DGbzenDmYqCQ5rT03qnG41INBhgkCQuhBCikCjySTytY5vBZEKS0QhjUpKDIxJCCCFyR5J42i1mRiNSjEaYkpMdHJEQQgiRO0U+iac3pxsMSDaZYNR7qwshhBAFXZFP4unN6WYzUsxmmCWJCyGEKCSKfBKHnsSV0YhUJyc42WleWSGEECKvFfkknj7sqskEi7MznPSkLoQQQhR0RT6JM0NN3OLsDBdJ4kIIIQoJSeIZeqdbXVzgIkMKCiGEKCQkiWfo2AYXF7hKEhdCCFFISBLP0JwONzeYAFDuFRdCCFEISBLP0LENbm4AgORbtxwYkRBCCJE7ksQzDLuq3N0BAAk3bzoyJCGEECJXJIlnSOIGPYknSU1cCCFEISBJPEPvdGOxYgCApMhIR4YkhBBC5IokcasVgNY7XZK4EEKIwkSSeIaauLl4cQBASlSUI0MSQgghcqXIJ3FkuCZu9vQEACRLEhdCCFEI5FkSV0rNVUqFKaUCslnfQSkVpZQ6pj/G51UsOcnYnO5UogQAIDU62hGhCCGEEHfFlIf7ngdgJoAFOZTZRfKpPIzhjjI1p+s1cUtsrCNDEkIIIXIlz2riJHcCiMir/duN3pxuNJvh4uWlLZIkLoQQohBw9DXx1kopP6XUOqVUA4dEoDenK5MJLiVLAgAoSVwIIUQhkJfN6XdyBEA1krFKqR4A/gFQ21ZBpdRgAIMBoGrVqnYNIuNgL66lSgEArHFxdj2GEEIIkRccVhMnGU0yVn++FoBZKVU6m7JzSDYn2bxMmTL2DSRjEi9eHBYAkCQuhBCiEHBYEldKlVdKKf15Sz2W/B+0PEPvdJPZjHgASEjI9zCEEEKIu5VnzelKqUUAOgAorZQKBTABgBkASP4IoA+AYUqpVAAJAF4g838y74zN6QCQoBQMiYn5HYYQQghx1/IsiZPsf4f1M6HdguZYWZJ4osEAJUlcCCFEIeDo3ukOl14TN5sBAMkGA4xJSY4MSQghhMiVIp/E066JG/UknmQySRIXQghRKEgST+vYpjenJ5tMMKWkODIiIYQQIlckiWdpTk81mWCWJC6EEKIQuGMSV5oq+RGMQ2RpTk9xcoJZH09dCCGEKMjumMT1277W5kMsjpFWEzcatZdOTnDWlz2ISCIgIABW/cuLEELYRXIykP93CRd5uW1OP6KUapGnkTiK1QoLAGXQToXF2fmBTuI//fQTGjVqhMmTJzs6FCHEgyAoCGjcGHB1BRo1ApYuve9d7vvqK1w+duz+YysCcnuf+CMAXlJKXQAQB0BBq6R751lk+cVqhRWAURs8DnRxgfMDWks9duwYRo0aBVdXV3z22Wd47LHH0K5dO5hMjhxCXwhRmEUNGQIVEICfDQb0PHsWtZ9/Hjh9GjAaAR8foHv3u9rfvjfeQOu5c3HUzQ1lIyJgdnbOm8DvQnR0NLZv344ePXoUuP+Xua2JdwPwEIBOAJ4G8JT+s/CzWJAxZdPVFa4PYJPQtWvX0KtXL5QqVQpHjhxBqVKl0KlTJ1SpUgV+fn4AAKvVitDQUAdHKoQoDCL27kXcnDnw3LkTP7i54eJbb6FfzZr4GwA++QT48EOgRw/ggw9yvc/g77/HI3Pn4oLZjCbx8VjXP8cxw/LFmjVrUKdOHfTq1Qsvv/wyUvSOzw4YYNQ2krl6AGgM4G390Ti329n70axZM9rT9pYtGae1KpAkN7dvTwK0pqTY9TiOlJqayjZt2tDNzY2HDh0iSV6+fJnz589nlSpVWLJkSa5fv559+/alwWDgvn37HByxEKIgSE1N5aZPP+XOli0Zf/Vq+vKwgweZpF0B5zWAm//5hySZmJjIDm3a8EWDgU08Pbn74YdJgNyy5Y7HsoSFMdxsZqDJxJuXLjGoVClGATw0ZUqevb87CQ4OZkV3d35VoQJP1KzJagC9vLxYq1YtFi9enEuWLLltG4vFwtTUVLvHAuAQbeVmWwtvKwSMBBAAYLL+OA7gf7nZ1t4PeyfxHc2bMyZjEn/iCRJgUkSEXY/jSLNmzSIALliw4LZ1ISEhrF27NgEQAN3d3dmyZUtaLBYHRCqEcKTk5OT0BJSUlMShjRszWk/WIe7ujD58mCS5r0ULJgH85cknOX/8+Ez7CAsL44gRI9inTx86A7zh7k42bUrm9D8lKYmXGjViEsCVkyaRJKMDA3nSxYUEGF2sGOM+/TRv3nQ2LBYLO7VsSX+DQUuVZjOjqlfn8Nde47PPPMNWPj50cnLi7t27M223dOlS1qtXjxcuXLBrPPebxP0BuGd47Q7APzfb2vth95p406aMypDEt/TuTQK8deaMXY/jKNeuXWOJEiXYuXNnWq1Wm2USEhL4xRdfcPbs2Zw/fz4B8OOPP5ZELkQREhMTw9q1a9PNzY1dunThK7168QrAW2XKcOfQoYwAmADwyCOPMA7gpipV7rjP9957jy/pXwL422+2C1mtjOjalQQ47aGHMv3fCQ0K4uyHHuIWfR/Jixfb6d3e2bKlS/kXQItS5L//kuvWkQYD6exMli5Nq8nEj0qW5COPPJK+TWpqKuvWrcv69evbvTZ+v0n8OACXDK9dABzPzbb2fti9Ju7jw8gMSXzbq6+SAK/t32/X4+Qni8XCdu3aceTIkezevTtdXFwYFBSU621ffvllAmCPHj1448aNPI5WCFEQjBgxgp0B/tW0Kd8vWZIz05LvgQMkyWPr13Nj2bKMAZgI8ICNlr2sYmJiWMrLiye8vMhSpciwsNvKJCxcSAL8wtOTVzM02Wf016JF3AswyWQi/fzu743mgsVi4bCqVUmAlqlT/1uxfTs5diw5aBDZqROtSrE3wKNHj5Ik582bRwBcvnSp3WO63yT+DgA/ABP1xzEAo3Kzrb0fdk/i3t68qVT6653Dh5MAz69bZ9fj5KfDhw+nN48D4MyZM+9qe6vVyh9++IFms5lVq1bl2rVr8yhSIURBEBAQwLoAU9OajvWHZcCA28qeDwnhrn//zfW+J06cyPoALSYT+eSTZGLifysTExlZsiT9AO7cti3H/bz36qu8BDC2TBlST5p5ZdmSJTwKMLpcOTI52Xah+HimeHszWCkOGzKEJFm/fn3+r2ZNWps3Jy9dsmtM95zEofVgbwOgKYAR+qPJnbbLq4fdk3jDhgzPkMT3vv8+CfB0Pjbb2NsXX3xBAPzss8/47rvvZtuMficHDhzgww8/TAAcPHgwk7P7MAshCrXx48dzDUBLsWLk5ctabferr0g79A26efMmvby8+E1aJ7fHHiN37iStViZ+9hkJ8P2mTe+4n8TERA5s1Ei7/AlwU8mS7PPkk/zggw947ty5+44zo4/r1tW+xCxcmHPBRYtIgM+5uHD37t10AxhVsiRZrx6ZlGTXmO63Jn40N+Xy42HvJL6zfn1ez5DE9+sfqoAff7TrcfJTx44d6e3tbZd9JSYmcty4cQTATp068fz583bZrxCi4HipVi0tHXz5ZZ7s/4cffiAA7h02jCxeXDvWwIFMcHHhGiDXd8RcvXqVr/XsyQW1atEC8IirKysbjfTw8OD06dN58eLF+4712LFj3AowsmRJ8k7XtZOSmFy6NNcArFy5MqeltWLs2nXfcWR1v0n8KwC9AajclM/Lh92TeN26vGowpL8+MnMmCfDotGl2PU5+iY2NpZOTE8eOHWvX/f7222/08PCgu7s7fX197bpvIYTjnD17lpMBWgwGMjIyT46RmprKFi1a0Gw28/+mTOHZXr1IgKkABz/66L3tdNky0tWVKWXKcFjz5gRAs9nMpbm8Hm21WhkUFJSpI53VauWobt1IgHETJuQujokTSYBNAcYbjeTLL9/Dm7mz+03iMQCsAJIBROuvo3Ozrb0f9k7iu+rU4eUMSTxA72Rx8OOP7Xqc/LJ8+XIC4IYNG+y+7/Pnz7N69eqsVq0aIx6gW/CEKMpmzJjBAwDj7fy/NauIiAh27do1va/Ou0YjhwPcfz+diI8fJ2vWJD08eHHePLZt25YGg4H/6PetZ+fMmTPs0KEDAfCNN96gxWKhxWLhmDFjOBtgitFosxOeTVeu0GI08kKWjoD2ll0Sz80sZgYAT5A0kHQiWZxkMZLF77RtoUDCqg+5CgBOJUoAAFKjox0U0L2zWCyYOHEiatWqhY4dO9p9/9WqVcPixYtx+fJljBkzxu77F0Lkvz3//INmAFyfeSZPj+Pl5YX169cjMDAQGzZswKH27ZEwcCBatmx57ztt2BDYtQuoUgVVBgzA1gYN0LJxYwwaNAjh4eE2N0lMTES/nj1R/cAB/FOnDtb++iu6deuGzp074+rXX2MoAOPQoUCZMrmLoUIFGPr2RVUAbNIEaN783t/PvbCV2bM+8ABfE99dsyYvGo3pry/s3k0C3G2jV2ZBN3fuXALIdXPSvRo5ciSNRiODg4Pz9DhCiLyVkJDA18xmrQZ58KCjw7l30dHkiBEkwJiWLVnKbGbLli05e/bsTB1y4+Li+Obrr3NDhh74Z2vVYo3q1dnVzY2pJhOtHTrcfae0vXu1/f38s53f2H9wrzVx3RalVG+lMlRZHxQkLBnelouXl7Y4Ls5REd2ThIQEjB8/Hi1btkSfPn3y9Fjjxo2DyWTC559/nqfHEULkrT179qBDSgqSixUDmjZ1dDj3rlgx4NtvgYUL4XHkCA43bIiw69cxbNgwPPPMM4iLi0NwcDDq16mDVvPm4XEAmDUL+P571AgOxtlHH8UGNzcYa9SAWrECcHK6u+O3bq3N5vbGG3nx7nKU2+lYhgAYDcCilErEf7OYFfomdWW1ghmTeMmSAABrIUvi3333HUJDQ/H7778jr79rVaxYEUOGDMEPP/yAt99+G02aNMnT4wkh7k5CQgLmz5+PDRs2YMCAAejZs6fNcps3bcJQAKpzZ8CQ2zpdAfbyy8D166g2dizOTZiAFRYL+n32GXr27AmnlBTMu3YNHQBtgpZhwwCrFdi3D/jnH6iyZYE1awA9B9y1hx+24xu5C7aq51kf0O4VfwXAeP11VQCP5GZbez/s3Zy+t0oVhpjN6a+TExNJgDs6drTrcfJSbGwsPT09+eSTT+bbMSMiIliuXDk2bdqUKQ/QZDFCFHbnzp2jj48PnQB6e3hQAfwxm1tmn27YUGsGvssBoQo0i4XU58AgwOu1a7MxwC/TmtB//dXREd4T3Gdz+g8AWgFImxcuBsDMnDZQSs1VSoUppQKyWa+UUt8ppYKVUv5KKYe05agsHdvMzs6IB4D4eEeEc0/WrFmDqKgojB07Nt+O6eXlhZkzZ+LIkSOYPXt2vh1XCJG9+Ph4dO7UCb1OnkS8yQS/2FhsqVABo0eMwIULFzKVDQgIQOkA/d9zHnSEdRiDAVi9GvDzA77/HmVjYnDQyQmjlQIHDQIGDnR0hHaV2yT+CMnhABIBgGQkgDtdNJgH4Ikc1ncHUFt/DAbgmEyQJYkDQKJSQGKiQ8K5F0uXLkX58uXRrl27OxcODwfefBM4e/a+j9u7d2907twZEydORGRk5H3vTwhxfyZNmoRXz53DxKQkGJ96Chg5Eh2vXsU0iwXjxo3LVHbGjBnobDTCWro0UK+egyLOI0Yj4O0NvP02cPQozD4+MJYvDzVtmqMjs7vcJvEUpZQR2v19UEqVgXbfeLZI7gQQkUORXgDSRtD3BVBCKVUhl/HYTdZr4oCWxA0JCfkdyj2JiYnBmjVr0KdPHxiNxjtv8L//Ab/8Atih1q6UwjfffINbt27hww8/vO/9CSG0v2mt9fTubNmyBYe++goTAGDAAGDFCmDGDGDoUAy3WuG/dCkOHToEAAgLC8MfCxeih4sLDB07Ag9gn+V05ctr171Pn773690FWG6T+HcA/gZQVin1GYDdAO63a3IlAJcyvA7Vl91GKTVYKXVIKXUou3v/7lXW5nQASDQaYUpKsutx8sqqVau0+x779btz4eXLgSVLtG/df/8NLF4MnD//33qLBQgJuavje3t7Y9SoUfjxxx+xePHiuwteCJFu7969aObtjWbFi6Ne7dr47rvvYLFYcrXtoUOH0L9nT8w3GmGpUQOYOfO/jmqTJ0N5eGCa2Zx+R8mUKVPQJSUFXnFxQO/eefWWCg6DAfDwcHQUecPWhXJbDwB1AQwH8DaAerncpjqAgGzWrQbQNsPrLQCa32mf9u7YdqBsWQa6umZaFuDqyoPlytn1OHmlc+fOrF69+p3n/j54kHRzI1u21CY1KFdO6+RhNJJps5S9+aa2bMuWu4ohOTmZbdu2pbu7O69cuXKP70SIoisiIoIjvLyYoBQJ8KqzM98A2OGxxxgdHZ3jtnFxcaxduzbnenjQqpQ2uUhW+tCgjQD+9NNPNJlMDKhcmaxYMftZukSBgvsZdvVeH3dI4j8B6J/h9SkAFe60T3sn8YNlyjDAzS3zshIl6F+8uF2PkxdCQkIIgJMnT8654JUrZPnyZPXq5LVr2rLgYHLNGtLHh/TwIF95Rfs4mM1knTpkfPxdxRIcHEyTycS33nrrHt+NEEXXB716MQ5gjLc3OXs2rW3akAD/VIrPPvVUjl/S33nnHbZL63k9cqTtQhERtLq58U8XFwJgczc3rfyUKXnzhoTdFcQk/iSAddDuOW8F4EBu9mnvJH6oVCn6u7tnWra7YkUGOzvb9Tj2lJyczOHDh7Ndu3ZUSuU8c09Kijb1n6sr6e9/+/rLl8n27UkXF7JDBy2xA6STE9miBfnRR9qUhFu3avvKwbBhw2gymXjmzJn7e4NCFCEXLlzgLoAxrq7k1avaQotFS7AApwH88MMPbW67Z88eOgG84uVF1qhBxsZmf6C336bVbObmBQsY3aQJWaJE7scHFw6X70kcwCIAVwGkQLve/QaAoQCG6usVtFvXQgAcz01TOvMgiR8uWZJ+Hh6Zlm2rU4fXM0yKUtBs3Lgxfcae559//vYCVivZpQv5ww/kd99pv+Z583LeqcWibUdqiXzcOLJVK6bda0lAm0Lw8cfJU6ds7uLKlSssVqwYW7duLXOPC5FL04YMIQFGfPDBbeus+rqnAf7555+Z1kVHR7NOnTr8rEQJ7e9z9eqcD3TmDGkwkJ6eWvlffrHjuxB5zSE18bx42DuJH/Hy4rEsTedbmjVjvDYiXYE0evRoOjs7Mza7b91nzzK9adzTU0voaQn6biUmatfQly8nhw0jS5XSmuZPnLBZfPHixQTA9957796OJ0QREhcXx0murtrf6/nztxdISKClaVPGGY1s7uTEg/r45ufPn2fDhg1Z1WBgiosL2bNn7g64dSv59NPa5bN7/Z8gHEKSeDaOlijBI56emZZt7tyZBJgSE2PXY9lL3bp1+fjjj2dfYMkS7Vfr5qZ987bVjH6vTpzQOsXVqZNt092QIUMIgOvWrbPfcYV4AE2YMIFHAEY1aJB9odBQppYvz0tGI73LlePKlStZvXp1lihRglcfe0y7FHb2bL7FLBwjuyT+AAyWe39s3SduKFUKABCV8farAuLcuXMICgpCjx49si904ADg7KxN0ff330CjRvYLoF497da0M2eAkSNtFvm///s/NGrUCK+88gpCQ0Ptd2whCrjU1FT8888/+O6773D58uUcy547dw7rPv8cTQAUHzQo+4KVKsG4ejUqms2YGxmJF3r1QomwMJz18UH5HTuAjz4CatSw7xsRhUaRT+IG8rYkbipdGgAQc+mSrU0cat26dQCA7t27Z1/o4EHAx0eblSibiQ/uS4cOwAcfAL/+CnzzzW2rXV1dsWzZMiQmJqJv375ITk62fwxCFEDDhg3Ds88+i5EjR6Jr166Ijo7OtN5iscDf3x/Xrl3Dc88+i6lWKyxeXsDrr+e842bNYPjjDzRNScHpKlVw0GiE19GjwIQJQJaR2ETRUuSTuLKRxJ3KlQMAJFy54oiQcrRu3TrUrFkTtWvXtl3AYgEOHwZatszbQD79FOjTBxgzRhunOIuHH34Yv/32G3x9fTFmzJi8jUWIAmDRokXY9ssv2NasGQ59+CHOnjqFN/SpKUliwYIFqFOnDho3bowKFSqgekAAOlosME6aBJQocecDPPcc1O+/o/LlyzCVLw8cOwZMnHj302aKB4utNvaC/LD3NfEAd3ceKF0607Jjc+aQAA/f6f7rfJaQkEBXV1cOHz48+0LHj2vXwxcuzI+AyCZNtM5uly/bLDJmzBgC4O+//5738Qhxn/7++2+2a9eOtWrW5LJly3K9XUpKCluWL8+rTk5Mu5vjRP36BMDdu3ezT58+BMBHmjfn0vHjOfW115hYvDjZoAGZlHR3QQYEkJGRd7eNKPQgHdtsO+Hmxv1lymRadka/V3rf22/b9Vj3a8OGDQTANWvWZF9IH5mJISH5E1RQkNaBrkoVcv7821anpKSwXbt29PT0ZJjckyoKsIMHD7KR2cztbm60APwR4ODnnuPJkyfvuO3yZcu4F2CSu7s2OuJnn5EAP3N1pYeHh9bRs1cvWitXTk/yLF4829s1hchKkng2Trq60rds2UzLrum12Z227sF2oFGjRtHZ2ZlxcXG2C0RFkV5eZK9e+RoXd+3SBoYByKlTb1sdGBhIo9Eoo7mJAishIYEPV6vGcyYTLSVKMLVfP6YqxQiAbwN8pX//HIc/fcfbmwRomTVLW2C1ki+8wFSjkXUA+tepo/19tGtH/vwzOX48uWNHPr078SCQJJ6NIBcX7itfPtOyhOhoEuD2Ll3seqx7FRERwa5du1Ipxaeeeir7gl98of1KDxzIv+DSpKaS/ftrx58587bVw4cPp9FozFWtRoj8NnfuXI5PqyFv3qwtDAxkUvv2JMCjACcMGWJz26CgIG4FGF28uDauQporV2h1c2NysWLafqdP1wZVEuIeSBLPxmkXF+6rUOG25TEAdzRvbtdj3av333+fSilOmDAh+ybpsDBtYJfu3fM1tkySk7VBJ2yMEBcWFkZ3d3fbI8wJ4UBWq5Xt6tVjglK09uuXdSW5bBmTDQbOVYp+fn63bT9r6FAS4K1PPrl95+PHa38PY8fmUfSiqJAkno0zzs7cU6nSbctDDQburlPHrse6F9evX6ebmxv79++fc8GhQ7UZyQID8yew7CQkaCPEGQzaKG8ZfPjhhwRAf3sOPiPEfdq8eTNHptXCs/lsxg0dSgvAFgCbNWvGvXv3ktS+APxcogRTlCLDw2/fMDmZXL9ea6kS4j5IEs9GsJMT91SufNvyU87O9K1Y0a7HuheTJk2iwWBgUFBQ9oWuXNGSZkHpiBcbS7Zpow37un9/+uKIiAh6enqyS5cutMqQj6KAeLJHD540Gmlp0SL7QjdvMlUfo9zfbGY3gFOnTuXWTZsYCvC8j0/+BSyKpOySeJG/T9xAgobbT0O8szOc4+MdEFFmvr6+aNiwIR5++OHsC/31F2C1AsOG5V9gOXF31+4dr1gReOEF4NYtAICXlxc+//xzbN68GUuWLHFsjEIAOH36NCLWrkVdiwWGwYOzL1iyJIxHjgBffokGlStjPQD399/H0q5dUQlAqREj8itkITKRJA7cNtgLACS6uMA1KSn/A8rCz88PjRs3zrnQsmVA/frao6Dw8gKWLAEuXdISeUoKAGDIkCFo0aIFRowYccdhKYXIazNmzMD7SsHq6al9TnNSowYwdiwMJ0+Co0fjbQCzAaTUqweP/v3zI1whblPkk7giARs18RR3d7jricdRwsPDceXKlZyT+LVrwM6dQN+++RdYbj3yCPDjj8CGDcCoUQAAo9GIBQsWICEhQYZkFQ517tw57P/5Z/QkYRg9GvDwyN2Gzs5Q33yjtYAtWwazvz/g4pK3wQqRjSKfxLNrTk8tXhzFU1MdENF//Pz8ACDnJL50qdYlpyAmcQB44w0tgc+aBezdCwCoW7cu5s6di3379mHs2LGOjU8UWePHj8f7ViusHh7AvTSHP/usNvSwyWT/4ITIJUnigM2aOD09URyA1YE1xRyT+MGDQHy8VtNt0QJo0CCfo7sLn34KVKoEvP22NrY7gL59+2L06NH4/vvvsWDBAgcHKIqa5cuX49Dvv6O31QrD//6nXf4RohCSJE7avCauKlQAAEQEBeV3SOn8/PxQsWJFlClTJvOKwEBtgpMWLYCTJwtOh7bseHgAX38NHD0K/Pxz+uJp06ahY8eOGDhwIFauXOnAAEVRER4ejmnTpuGll17CN2XKQLm6AqNHOzosIe6ZJHHAZk3cpVYtAMCNo0fzN6AMjhw5YrsWvnat9vPECW32o+efz9e47km/fkDHjtrcxzdvAgDMZjNWrlyJZs2aoU+fPpgzZ46DgxQPsoSEBLRq1Qrvv/8+BjdqhCciIqCGDgWyfkkWohCRJJ7NNfHiek/v6BMn8jskAMCOHTsQGBiIxx9//PaVGzYADRsCv/8OzJ0LuLnlf4B3Syng+++B6GitWZ0EABQrVgwbN25Ely5dMGTIEPTr1w/Xrl1zcLDiQTR9+nRcOXsW+6dOxfdXr0JVqwZ88omjwxLivkgSB2zWxMs0aQIASAoJyd+AoA3A89FHH6FixYoYMmRI5pVxccCuXUC3bsBLL2mdawqLBg2AyZOBxYu1a/k6T09P/Pvvv5gyZQpWrVqFNm3aIDQ01IGBigdBamoqrly5gk2bNuGdd97BrM8/R3Dx4mj5/vval8m//5Zr4aLQK/LdKg2AzZp46Tp1kADAeulSvsbj7++P4cOHY8+ePZg9ezZcXV0zF9i+HUhO1pJ4YfTee8Du3Vpt3M0NeO01AIDJZMJHH32Exx9/HF26dEGXLl2wf/9+eHp6OjhgUVhERERg1apVWL58OQ7v34+HbtyAAUA9ADUNBmwuVQoVIyO1Fqzu3YGSJR0dshD3rcjXxLO7T9xgNOKayQSn69fzNZ733nsPJ06cwIwZMzDY1ghSS5YAxYoB7drla1x2YzBot8V17gy8/jrw3XeZVrdo0QL//vsvgoODMWzYMG1sYCFyYLVa8fnnn6N8+fIYPGAAmu7di4DEROwGsBPATwDGAWh08ybUjBlaC5YkcPGAyNMkrpR6Qil1SikVrJR638b615VS4UqpY/pjUF7GY4sRsJnEASDS3R3ukZH5FktycjJ27tyJF198ESNHjoTBYAAuXwbGjweeeQbYsUNLgK+8UrgHl3B3B/79V7sUMHIkMHiwNrKbrn379pg0aRIWLVqEH374wYGBioKMJBb+9hta+vhg50cfYWfZsogrWRKTIyNRqkkT7Qvvxo3AqVNQycna8L/Dhzs6bCHsy9aA6vZ4QMuPIQBqAnAC4AegfpYyrwOYeTf7tfcEKHEAt2Uz8cHO6tV5yWSy6/FysmvXLgLgX3/99d/Ct94ilSLd3LRZynKYaanQSUkhR4/WJkoxm7WZ2C5cIEmmpqby6aefplKKK1ascHCgoiBa9MILvJE2+xhAa4UK5IsvkmvXalOICvEAgQMmQGkJIJjkWZLJABYD6JWHx7sn2V0TB4CU8uVRLjUV1AcoyWtbtmyBUgodOnT4b+HRo0D79tp1ZGdnoG1boFGjfIknz5lMwDffAMHBwKBBwK+/ArVqAZ98AqPBgEWLFqFly5bo27cvvvnmG1itVkdHLAqA1IQEbGveHC8sXozoUqVgnTwZWLgQ6vx54I8/tOvdNsZ+EOJBlJdJvBKAjL3CQvVlWfVWSvkrpZYrparY2pFSarBS6pBS6lB4eLhdg8ypOd1QpQrMAG7m021mW7duRdOmTeGV1mPWYgH8/QEfH6BJE+D4cW285gdN1arasKwhIdokFFOmAAMGwN3FBZs3b0avXr0wZswYtG7dGr6+vo6OVuQBkoiMjERUVFS2ZaKiovDT0KE44eWFjocPY1ejRqh68SIMn3wCvPwy4OSUjxELUTA4umPbvwCqk/QGsAnAfFuFSM4h2Zxk89tGL7tP2d1iBgDODz0EALhx7Jhdj2nLtWvXsHfvXnTp0uW/hSEh2i1laQO+1Kz5YA9MUaUKMH8+MHGi9nPsWHh4eGD58uVYsGABLl26hNatW6NPnz7YunWrdHp7AISHh2PUqFGoXLkySpYsiRIlSuCZZ57B2rVrER0dDUDruDZ/3jxMq1wZb/z0EyqT2Dd6NNr5+8NYGMZIECIP5eUtZpcBZKxZV9aXpSN5M8PLXwBMz8N4bMqpJu6pj0eelwO+fPnll/D390fNmjWRmpqKgQMH/rcy7cuDj0+eHb/AUQqYMAGIjARmzACqVYNh1Ci88sorePbZZzF16lTMmjULK1aswNNPP405c+agfPnyjo5a3ANfX190794dsbGx6NWrF9q0aYPw8HD8+OOPWLlyJVxcXPDEE08g6ORJDDx1Cp8DuNW+PUquXInWJUo4OnwhCgZbF8rt8YD2BeEsgBr4r2NbgyxlKmR4/iwA3zvt154d26ypqSTArY89ZnN9+KlTJMAdTz9tt2Nm5ePjQwA0GAzs3r07mZRERkVpKz/4gDSZyMTEPDt+gZWaSvburXXqW7Ik06qEhAR+9dVXdHV1ZZs2bZiamuqgIMW9Cg8PZ+XKlVmzRg1e+vprcsQIrRPn3LmMDwvjli1b+NZbb7FihQr8s1o1rePasGGkxeLo0IVwCGTTsS3Pkrh2TPQAcBpaL/WP9GWTAfTUn38BIFBP8NsA1L3TPu2ZxFMTE7Uk3qmTzfVWq5XXleKuOnXsdsyMUlJS6OzsTA8PDwLgunXryMGDSS8v8uRJsnt30ts7T45dKMTHk23bkgYD+X//p33ByWDhwoUEwGnTpjkmvgdQSEgIP/jgA3bo0IFPPfUUp06dyqi0L5V2YrVa+cwzz9DNyYlhfftq/4bc3clixbTnpUuTX35JHj9O9u+vLXvzTUngokhzSBLPi4c9k3hSTIyWxLt0ybbMkeLF6V+smN2OmVFQUBAB8Oeff+aBAwe0W668vLRfi5eXlrzeeCNPjl1oxMaSPXtq56RkSe1LzvbtpMVCq9XK3r1702g0cvny5Y6OtFCLi4vjpPfe4ztGIw8CTDAY6Ofqyg8BtvLy4vHjx2/bxs/Pj7/88gv/+ecfJicn5/pYK1asoDvA4Hr1tN/ruHFagrZYyD17yG7dmHbbGA0GcsoUuWVMFHmSxG1IiIggAW7r1i3bMtvr1WOkUrTasRZgtVqZmJjI5cuXcxLAm48/rh9s+3//1Hx8tJ83btjtuIWWxUKuXq3dA+zurp2jKlXI2bMZExHBNm3a0Gw2c9q0aUwsipce7oLVauXfixfz81at6LtlC61WK9f+9BN/KVaMN/XEmdSkCTl8ONmyJQkwFeCAUqW4atUqBgYGcuvWrezRowcBpD+aNWvGFStW8NSpU0xJSbF53KCgIE6bNo3lvLx40MODVoOBnD3bdqB79mjrLl3K4zMiROEgSdyGuLAwLYn36JFtmW3PPksCDLdRE7lXy5cvp4eHBwcMGED/tBpHcDD5zjukkxMZE2O3Yz1wYmPJRYvINm2089a3L29FRvLZZ58lACqlWLFiRbZu3ZrPP/88e/fuzY4dO/K5557jjBkzeOXKFUe/g7uWkpLC2NhYxsfH848//qCvr+897efy5ct8tXVrHtI/c5cB7nRyYgLAZKV4vVMncseOrBsxtkEDxgJ8OkPS7lm8OA8+8ghjO3ZkUPfubO7llb7OycmJDRs25Ntvv83Vq1dzypQprFmzJgGwIsB15ctrv7uFC+1wdoQoGiSJ2xB9+bKWxJ96KtsyB7/4ggR49Ouv7XbcwYMHEwDNAJPTkvi4cWT16tp1cHFnViup/2744Ye0Wixct24dJ0yYwAEDBrBTp0586KGH+PDDD7NNmzasVasWAdDV1ZUff/wxo6Oj8zXcJUuWsH///gwLC8tV+ZSUFC5YsICPP/44S7i6sixAN2dnVtA/N4MGDeLp06eZkpLCQ4cO8bs+fbjGZOJcJydO6taNSVn6D6z65x+OcHdnDMAEd3fGfvEFT1SpwjOlS/NE585MOXMm+2CuXaPl4YdJgBHVqvFaw4baeXdzIxs1Il1caHVy4rW+fblq0iS+O3Yse/ToQWdn5/TE/ra3N883asT0ZvL337+f0ylEkSNJ3IZb589rSbxXr2zLXD5wgAS4o18/ux23mY8PqwOsn/YPzclJ+6kUuWaN3Y7zwLNayQEDtHPXpQt54AB59iy5ciU5c6bWOervv7UOciRPnjzJF154gQBYrlw5rl69Ol/C9PX1pZOTEwGwRo0aXLp0qc0m5zShoaFs27YtSwP8vVgxxumfD4v+eblRvDhfNBpZGuCrAHfoyyOdnRlnNjMR4EdNmnDBggU8ePAgv+zbl8f0MnEtWtxbE3VSEjlrFtm+PVmpEjl9OpmQoK27fJkcOPC/z3HZsuQLLzBmxgweHz+eyU2bMr3D2vjxZEDAPZ5JUZil7txJy9Gjjg7jvpz4+WdGnDrlkGNLErchMiRES+LPPpttGavFwiiAOxo1sssx4+PjOchgYIpSfD8tiX/8sfYzu+uDIntWK/nTT2Tx4v/V8rI+nJ2167s//kgmJdHX15c+Pj40GAz85JNPGJAHScVqtXLVqlXs2rUrzWYzXy5fnmd792a96tUJgC1atLDZtL9x40aWKVOGPVxcmODhQavJRL76Kvntt+RHH2m99OvXz/T+osqVY8z48WRcHBkRwbAqVbRkD/Cg3toTUawYk+fPz9se3uHh5Ny55EsvkWlN5gBZowb5ww9afKJASg4PZ8ATTzDp2rXcbWC1MmzcOCYeO5a74hERjDYaeaVYsdx/Bq3W9C/geSFs3z6GFC/OMz//nKvygV9/TQtAP09Pu/aRyi1J4jbcCAoiAW7v0yfHcv7u7jxWvLhdjrl3717+pv9zS/TwoNVo1O4DDw+3y/6LrKgo7UvQrFnkvn3k1avkrVvkxo3kmDFks2bax71aNfLnnxmb4To6AP7000/Z7jo2Npa//vorN2zYkKtQrly5wo4dOxIAq1atyhmvvEKLqysJ0NqmDVdPnUp3d3dWqlSJa9eu1cOP4pAhQwiAIypXpsVsJhs00G6zyiolhdy0ifz0U+0adtae2xERTJ0+nRH9+vFao0a89vTTZEREbs+kfVitWo1740YtXpFZQoI2FkIeiQ4KYsCYMbnu1b//iSdIgL7PPJOr8pe++44EePyhh25faeN9nXnppfQvddd//fW/FadPa3+zWZO1xcJzTZsy3MOD1tz2EbJaGX0XX8iP6OMPBLu7p8ecevkyz3h7M7hOHSbu2UMGBpIBAbw8dSojDAZG6u/h6Pjx/x3WYtHGHMljksRtCDt+XEvizz+fY7ltLVowEWDs9ev3fcxvv/2WARlrifXq3fc+RS5YreT69ek9rlmjBrl2LUNDQ9mjRw8qpTh8+HAuWrSIf/31F6Oiorhu3To+99xz9PT0ZHmAngAHDBjA8+fP2zzErVu3OH36dJYsWZKurq5cOGUKU998U7t2XLWq9iXDxUWrPbdqxaf1zl5PPvkka9WqRYPBwCVdu2q9tlu3Jm/ezOeTJPKDNSWF511d6ZvN7In2sE9vrTmXcUbEbCSFhzNSKRLgRWfnTDXlmC1bGH/gQOYNUlIY6uHBtEs8CSdOaMtv3eLNVq0YXaoUrWlNzjdu0DpvHmOMRq53cuJFgKHlytH6zjtM7NqVqQYDCTCmZk1yxgzy99/JVat4tU+f9P+RAS+/nKv37Nu9Oy0Az2TpMGm9deu2LxYhU6ZoX1r0v8eTDRsypHp1xhoMTAAYYaNF77RS3PfbbzxlNjMZYJjJxJNeXrxqNDJKKe596ile/Osvxl+8mCe3REoSt+Ha0aNaEn/xxRzLHfr8cxLggU8/ve9jvt67t3Zts2xZpvWuFvnIatVuV2vQQDv/Xbsy9Zln+E7Pnpk6Yjk5OdEL4LuenjxctSqtBgMTnJ35kcFAb4OBtWrWZOfOnfn000+zSpUqrF27Nt3c3AiAT3TrxusjR2rXiJ2dtXv9z57Vjh8WRk6bRnp50Wo08lDr1qxevDiblCvH8PbttZgef1zrhS8eSMf0BBKhFC1ZOiBmZ3fPntzdtGmuyiaEhTFGTzx7Hn00fXnSlSs81rIlQ779Nn2ZNTGRRx55hAT4T/XqJMDz777LuL/+YujUqUwGmKAUwz79lPz3X1oWLuQlvZPjzHr1mALwgrc3Y157jdElSjAZ4E2AsS4ujPX0TE+AgQDnf/wxf9E7RcYD9Af4HcAhxYvzmo2kucTdnVvd3XnLaGTY1KmMWLCAh3v25EUnJx4aMSLTew47cIDx+nZHypdPX3525kxGK8VzxYrx/CefMKh3b57W3+cxk4lXz57lQScnxgHcD3CJpyf/mjyZO5Yt48+PPMIfO3Tgr088wSUjRvCqfvnLb9ky/luvHtdXrkxfd3duL1+ee0uWzBR7LMAbJ0/m6veVW5LEbbhy8CAJcMcdvunF37jBRIBb7+PYVquVCxcu5GNpv+iZM7WfU6bc8z7FfYiP14b69PYmPT3J0qWZ+NVXvDR5Mi8/8wxXtmnD6LR/QuXKaXcPdO+e/kca5ezMjaVK8aMKFTixUye+2Ls3hw4dyqMbNpD6bYns14+8eNH28cPDtWFGDQZavbxo9fTUEv7nn5N3MXCKcDyrxUK/zz9n3NWrmZbvmzGDm0aOvK38oQoV0jspnvjhh/TlZ7/5hicffpgpWS59hMyfn17+YoaOrxGBgfRt2ZI3/f0zld8/eLCWjA0GXjGZaE1NZdiKFQzVa52pAAMaN+bJatUYpndG/KdqVYZfvJg+VkDa46DZzP0mU6Zl4QC/qVmT0VFR/FfvixILcDXAd5o148Q+fbgK4G8A3zebOdzbm/N++40Wi4X+fn58q08fvjtmDGfMmMFDhw4xJiaG4z/+mMP69eOIbt04pn17Tnj1VZ48cYIbvvoq05zxBHhBbzXwq1ePRx59lIfbt+dFZ2fGAlxfqxYJ8GjjxjxTujQtAANMJl7St00BGADw94cfZmhICEny0oULPHzgAOPv4/q71WrlkT/+4NYRI7jxySe5qXFjptr571iSuA2X9+0jAe587bU7lj3i5cUTLi73dJx//vmHFSpUIAD++NBD2mkPD9cGd7l16572Kezo9Gky7bYpQGv+Bsjatcm9ezM3jZ09q3XeeuUVLbmnbVOtGvnyy9pIe2Yz+dVXuWtS8/fXRih74gktDlHwWK28vnp1ejNz2JEj9C9Thqd//50keXjCBBLg3gzXh4M3bOAtPWn4/fJL+vIgvXPUev3e+7167frGypXpNcmjbduml4/as4ehzs68YDAwHuA+fRhma0ICT+qjO57w8qI1NpaWwECeHjeO10wmnjGZuPG110iAl/Sm70sA/377ba719GQYwINKcV3x4lwxcGD6/ANrp0/ngldf5cJBg/hbly48degQg/z9Oeull/j9a6/xu2HDuGrxYlr1z/b+vXs5ddw4zvz+ex44cIAWfSTFAwcOMCQkhBY7dADbv3cvl37xBZe99x63/vwzI69e5dqyZXler/UnAPR1dubGESMYfvYsQwwGXgfoazZzSf36vBQUxNNHj/Kfzz/nwd27GVNIx+GQJG7DpV27tCQ+cOAdy+7Qh4IM3rTpro5x4cIFenp60tvbm7NmzmRK587a/eCiYLFatQR9/Lj2/PLl28Zqv43FQp46Ra5YQbZood1C1a8fmXaNUNyz4GXLuPcuLjXdOn+eO196iUm5uP8/9to1BukJOI01w3Yx/v488fjjTNKbTw/qyfCAfpvpfv2L+Gk3N1ri43nGxYWpegL2/+ILHv3iC55ydmakUrxiNPKSycRjlSrxot65Mchk4sWDB7mjTBkmAYzUa8PnleLq0qVJgCHVqvF86dJMAngd4KZJk7ipenXGAwyqXp3X9Vr1v/pdAGk1dQL0M5m4+5tvGH7+PPeYTNzm7MxFzZvzrN7pKyUlhVevXk1PxIVVcnIyb926ddt7iY6OLrSJOieSxG24sG0bCXDXm2/esWz4sWNMALilWrVc7z8qKoqtW7dmTxcXxrVrR3bsqJ3yCRPuOWYhHhSRJ09yf7t2TLLRc95Pv5RxccuWTMstV6/abOHYqzej7uzdO8djRgQHM1BvaTmoDzhzpGtXJgM8PX06rfHxDNGPHVijBmNPnGC0UkzRm4z9X3xRS+j6Ps7riXnjSy/xqt5JiwDDlOKhKVPo+/nnvALwuNnMbaVLc1XHjrylD/izc+ZMLi9dmsu8vDi3SRPuWb6cFwIDucLNjTsAbjSbubhSJR7R74o4sngxfQ0GHlSKazw8uPCpp5iamsoljz3GhdWq8dcOHbhqwgTGZPhCUtgTtfiPJHEbzm/aRALcPXRorsrvad2aqQB9MzSP2XLjxg1+99139PHxYQ2jkUnu7mSpUtq9s9Ony2QOhcCxsmW5J4dBgO5WwJw53NG5833v58Bjj/HYxx/f1z7OL1jAWyYTQ1etumPZ42PH8oyXF2MvXLivY9qyR+9cuPf11zMtP7lgQXoy3Jk2rwDJkDlzmArwSIZlJHnq22+Zdl02zGBgYlgYU0NDGbZ+PYO+/poHBw7k3qZNeVIfsS4R4BmzmTEAj+m16pv68it60/OqDJ2y4gGuHD06vbPYcZOJ4RcucEPZsjxmMnFVy5ZMTUmh3z//cFP//tw1bhzjIyPT47vbRCqJV9giSdyGkLVrtSQ+fHiuysecO8eb+vWW1R9+mGmd1WrlyZMnGRYWxgYNGhAAO5Qty6hq1bQpFnMa1rIQunL5Mrdt3HhX21itVs728eEGOw5ha0v4yZM8UqwYT69ceU/bX9JbaE66ueWq/JEpUxgZHJzteqvFwpN6jS10716bZSICA3mgf39ac7in+tScOSTAUznElbBrF+N37sx2vSUmhpf0pljfOwxgdGrKlPRmYt9c3uaTVdyVK9zXuDHPZrnV6crevelDDgd6eKQvt6am0rdSJd4CGODqymC9H0pMQACvG41M1puOg0eO5Jn//Y8BLVsyBeAZo5E7Jk5MT7xZH9EAD3p6cru3N/1nz+a53bt52NmZFwCuqViRx7dt42p3d651cuKvjz3G+Lg4/la/PufVrMkV+t/69lmzuH7aNEbm9z33QlCSuE3B//5LAtyT5XaFnET6+jLU1ZVRANd99VX68mkTJ/IlgJWU4o8GAy1Go3Z6S5fWBrwoQAICAu772/7yRx7hOYAxdzHX9IVjx0iAWypVuq9jpzn+yy/agD1Z7Bw4kAS4NUMHobuxS79HNRVgrH5dNHDCBAbPmXNb2YAfftBqjDnM+35k2rT/apbZJMPDlSqRAI/mcLeCn369lADDst67S/LWjh2MV4qRRiNTbIzRbrl2jSfr1NGOYzIxVikm27oX3WLhibffZirAQ87OPGE2MygX0/HGBQczfMWK9JYma0oKj1aoQAK8ZDYz4fp1MiGBV1ev5nFPTyYC3NSkifYl4a23uGfAAB7Ub73c0q4dNz/3HAkwqEwZJus14vXjxzM4Q7N1LMC/ypfnqX37aLVa+VfnzvynYUOu7NaN6wYP5q4vv+SJzZsZm821cqvVKjVfUShIErfhzN9/kwD3jh59V9vFnz7NG2YzzwP8tFYtDu3cmWuzfvt//XVy6lQyt8MYZiM6Opr+WW4huROLxcInHnuMC+fPv23dgf372Rfg8j/+uG3diRMnOGjQoNsmz7Blr97suDeHkc6y2vbhhyTAEIPhvnut3jxzhkkAd9aqddu6PfpITMeySTwRoaGZXidkGf70sN6hiACPTZ3KyMBAJgK8bDJlqilbU1N5Uj8P4QYDU21Mg5qakMBjnp68ZjAw1Gikb4UKt5U58eWX6Z+bgzbWk6Tfu++SANfpzb/7Xngh0/rIlSt51dmZ1/Saqn+nTmRgIJNWr+a1qVN5+oknGG00MhHgvCZNuHPqVBJgQN26DOjcmQGPPsrjzZvzeKNGvFisGAlwe7FivHLmDNfpo3mFfPABQ2fM4IXp03n6vfd4unRpXnN25qEuXbi/bVvG6Lf+HK9YkYebNk2/pWldzZraFyJ9PQFGAtzy6qsMP3mSCRn+buIAbnn2WVotFl45fpz7nZ15wGzmyoce4qGlS0mSR3bu5N/jx3Pb8uUMOnlSkrAoEiSJ23Bq2TItEY0de9fbRm/dyhv6LR5pj5RPP9XmYc7S8/V+fNmzJ6cD/CnD/aTh4eF88cUXsx3ze+/27bwM8CcbveDn6h1z5mYYBCLNt9278yzAzcuW5RhTTHR0+ohGq/RZ1ywWC29kM/f5jRs3GBoayjVp04cCPL1rV47HuJPtr7+engwSM9SyLCkpDFeKFmhjhkdnSdCbevViLMDAxYtJkvtefZUpAI9/8w1JMv76dSYB3FynDpMB7mrfnvsyxO03ebK2I6uVR556igS4Q689Hp06NdOxok+d4v6qVUmA215+mdvr1mUUwOTYWFojIxm1bx/9hwxhrFI8ZTJxg7c3UwGemjGDZ2bP5skvvmDA2LE83K0bkwHud3dndHg4T5nNPOvmxhPPPMMTbdrwkj7QxAWAayZN4r8Zb33TH0kA/3Fz41+ffkqr1cqU5GTu1e+tvQnwKsDzAE8D9DWZuLRnT0bpzcbnMgykkfFxTinu0ud3TwG4w8uLK1u14nV9f7vc3bnqxRdptVq54oUXuKZSJf7TvDn/7dOH5w8dSj9P/vPm8cC0aTy3bRuTZHx1IWySJG5D0OLFJEDfe50W0Wol9+/X7hvO0ovWHiwWC9c7O5MAvwY4ceJEJicns1+7dtwOcHyWDj5pZuvNwfsBRkZGMi7DP8Z5ZcqQAFe4udFqtabXiFNTU/m3fqy53brZ3G9kZCQ/++wzbpw9O/0f+cayZUmSC157jf8YjbykD6CQxmq18o8KFbjG1ZU79CZUAtyUy34I2dlfqlT6NdVNPXpwd6VKPL1yJU8tWUIC3FWjBglwT6VK9Hd355k1a3jVz4+39G0OeHnxyoEDjE5LSE5OTLh4kUf1pLv///6P/q6uvOTkxFiAW8qX5yWDgUHFi9O/b1+e0hPl+vLleevqVUYCPO3uzoMNGvBwzZo8pyc3Atykf9HZO24cCaSfg7SHr5sbj69bx5AtW9KvQWd8pALc4+rKq/pQlqv1LxXJ0OYE32Y0cq63N08cPkySDNy9mz83bsyfO3fm/Dfe4Irp07lz48b0e4HT3LhxgydPnuT58+d5/fp1RkdHZzu72s4//uCqTz/lmmnTuP7//o+bZ89muH5rT4i/P8OvX0+vEcfHx+eqNUcIkXuSxG04sXCh9k80Sye1/LT9iy+49n//s7nu0L59jAaYqN/OMgmg2WDgRv2f+0KDgZGRkbx48SJ9GjfmnDlzaLVa+YveszYV4IBnn6Wb2czly5czNDQ0fUrKiwAHDhjAihUr8uzZs9y+cWN6glvm7s41a9Zw3rx5meKZ9tJLDAL4jR7P1WLFeE4pRkRE8KDeB+CPLl1IkiHBwZwxeTLX/fRTerJNBHikTh0mAlzbqBGtViv3btnCwKNHeSM8nMsmT+ZJPz9uW7qUvw8ezMP793PRiBFc8tZbPLB5M5f17cvV48Zx148/MgngrqZNGZWhifaqwcBDpUppSXnz5vRkmQIwXKn0MY93N27MtFp8HMCNaSOsQWuKXvXEE7RarVzXogUJ8IjJxMDVq7mmU6f0cgFKcUmrVkzUp+Nc06KFVptViqeMRu709OSqVq3ot2hR+vmLj4jgRm9vrvf25uoOHbi6f39umzaN8RmGWPWdOZNb33uPOz/7jHu+/577//yTIYGBmZJrbEwM965bx+N+frx8+fJtyVkI8eCRJG5D4Lx5JMAD+Xzf9ok//+SuKlUYGRrKEKORMQDPHTlCkpw0diyfa9aMJ06c4C/6IBNRv/5Kq958nKwPgZhYoQJvAfxk3Dh2bdiQpwAOMxj43nvvMQhgnH5dczbAKIB9PDw4Rm9Kj9UHiJgKbfzi7k2a8MNWrbR1xYrxKsDXAI4F+MnHH3P06NH84osv+G+GhJmoFI+/8goJ8JNGjbTYlGKQUvzxu++4yNmZKQAP6DXGCD3uwwMHMtDLi2eU4s/lyzNGb3o9ou/7CrTOSgTSv1SkfSHJWkM9s3gx9zRuzEiluP+dd3hDKd6C1imKJLe1asXtLVvywvr19C9ZkgElStD3jTeYEhdH34ce4p66dXlsxgwtYffqxdVt23LH55+n/54irl3jniVL0muViQkJ9P37bx7fv99mTVOuzQoh8ookcRuO//KLlsQnTSKp3XK2t2JFnsvj3uQH9NrirgzX1Jc2bco/f/2Vh6H1wu1pNPILk4kpgDY0q9WqNdu/8w45dy6tq1eTAJ8EuErfR5xS7Jv2fPJkJqb1kIc23vDH+nOr3gKR9vgX4J8Ak41GRn72WaZ13wP8WynO01/H16xJAgwuXZrx+rC1KQBjjUZe0yd2iNLL3tDntQ5u25bnR4wgAUZs3szzH37IBD22kIYNGVyrFq94evLc0KE806ABT7Vpw7BvvuGZFi146aOPGDp9OgO7dWPE33/z/IwZDBwzhlf0e5wtiYlM0e/JvXXhAhMy3J8rhBAPCkniNvj/+CMJcE/VqtxVqxav6beu7LbzsKiX16/nwX79mBITw8DffyeB9IEjIpXi6YoVeQ3gBn1Zco0atCjFRKOREfXr295pUhItnp5MSRvne/RoWvTrsNbixclLl5jcuTOtSpFff639BLQJP1JTmVqqFC0lS5J6j3EC2mxbV67QajIxtW1bpurX1i1lytBiNjPZy4sMD2fiQw8xdswYLY516xjdsSOjp0whU1MZ98EHjOzfn0nz52vDkq5YoU2pmZqqjUOeJiVFG9pUCCHEHTkkiQN4AsApAMEA3rex3hnAEn39fgDV77RPeyZxP/0eX+pNuBdMJu6tVIkp+G9QDmvWW6GsVp6eN48HR43i2aVL0yeD9/vsM4a4uvLI6NHc364dT5YowZNffcWATz/lLT2BBru68qzZzFsAg77/ngTo26oVo/79lwnOzkx2dmb8xIlkVJQ2NOtTT5E5zQc8Zw75/PParWypqeSqVeRnn2nTXZLk0aOk3gubS5aQS5f+Nx74pk3kgQNaop07l9y377+R5E6c0Gb5SkoiN28mExO1RJzW09sOkxoIIYTIveySuNLW2Z9SygjgNICuAEIBHATQn+SJDGXeAuBNcqhS6gUAz5J8Pqf9Nm/enIcOHbJLjH7ffYfGI0ciGUBUYCBK162LK/v3o1ybNgh2ccGtEiXgc+0abppMuFK9OpKqVUMZX188HBeXvo/LZjPCS5ZEw+vXkQTAXV8eBqCs/vys0Ygzffqg9vLliHF1RfzQoWj95ZcIW7YMJTt3hqlkSbu8HyGEEA8mpdRhks2zLjfl4TFbAggmeVYPYDGAXgBOZCjTC8BE/flyADOVUop59c0ii+ToaABAQIkSaFq/PgCgUuvW2Pniiyj799+oGh6OA/XrA5GReDgkBOWCgxHg6orNvXujXJ8+CNuwAS7//ovit25hV+XKqLthA3YMHw5nb280//BDbPzkEziVLYsGb7yBbjVq3Hb8sn375sfbFEII8YDKyyReCcClDK9DATySXRmSqUqpKAClANzIWEgpNRjAYACoWrWq3QKM3boVABDdsWOm5e3/+CP9eUX9Z0pSEq6fO4eGdeuiYdrKF164bZ8Vtm1Lf/74nDl2i1UIIYTIyuDoAHKD5BySzUk2L1OmjN3222r5cuz/+GO0zEWyNTs7o1zdunY7thBCCHG/8rImfhlAlQyvK+vLbJUJVUqZAHgCuJmHMWXiWrIkHvn00/w6nBBCCGFXeVkTPwigtlKqhlLKCcALAFZlKbMKwGv68z4AtubX9XAhhBCisMuzmrh+jfttABsAGAHMJRmolJoMrav8KgC/AliolAoGEAEt0QshhBAiF/KyOR0k1wJYm2XZ+AzPEwFIF20hhBDiHhSKjm1CCCGEuJ0kcSGEEKKQyrMR2/KKUiocwAU77rI0styXXsTJ+chMzkdmcj4yk/ORmZyPzOx5PqqRvO0e60KXxO1NKXXI1lB2RZWcj8zkfGQm5yMzOR+ZyfnILD/OhzSnCyGEEIWUJHEhhBCikJIkDsgA55nJ+chMzkdmcj4yk/ORmZyPzPL8fBT5a+JCCCFEYSU1cSGEEKKQKtJJXCn1hFLqlFIqWCn1vqPjcQSl1Hml1HGl1DGl1CF9WUml1Cal1Bn9p5ej48wrSqm5SqkwpVRAhmU237/SfKd/XvyVUk0dF3neyOZ8TFRKXdY/I8eUUj0yrPtAPx+nlFLdHBN13lBKVVFKbVNKnVBKBSqlRurLi+TnI4fzUVQ/Hy5KqQNKKT/9fEzSl9dQSu3X3/cSfe4QKKWc9dfB+vrqdgmEZJF8QBvPPQRATQBOAPwA1Hd0XA44D+cBlM6ybDqA9/Xn7wOY5ug48/D9twfQFEDAnd4/gB4A1gFQAFoB2O/o+PPpfEwEMNZG2fr6340zgBr635PR0e/BjueiAoCm+vNiAE7r77lIfj5yOB9F9fOhAHjoz80A9uu/96UAXtCX/whgmP78LQA/6s9fALDEHnEU5Zp4SwDBJM+STAawGEAvB8dUUPQCMF9/Ph/AM44LJW+R3Alt8p2Msnv/vQAsoMYXQAmlVIV8CTSfZHM+stMLwGKSSSTPAQiG9nf1QCB5leQR/XkMgJMAKqGIfj5yOB/ZedA/HyQZq7806w8C6ARgub486+cj7XOzHEBnpZS63ziKchKvBOBShtehyPkD+aAigI1KqcNKqcH6snIkr+rPrwEo55jQHCa791+UPzNv603EczNcXiky50Nv+mwCrbZV5D8fWc4HUEQ/H0opo1LqGIAwAJugtTbcIpmqF8n4ntPPh74+CkCp+42hKCdxoWlLsimA7gCGK6XaZ1xJre2nyN7CUNTfv242gIcA+AC4CuBrh0aTz5RSHgBWABhFMjrjuqL4+bBxPors54OkhaQPgMrQWhnq5ncMRTmJXwZQJcPryvqyIoXkZf1nGIC/oX0Qr6c1A+o/wxwXoUNk9/6L5GeG5HX9n5UVwM/4r0n0gT8fSikztIT1B8m/9MVF9vNh63wU5c9HGpK3AGwD0BraZZS0ab4zvuf086Gv9wRw836PXZST+EEAtfWehE7QOhqscnBM+Uop5a6UKpb2HMDjAAKgnYfX9GKvAVjpmAgdJrv3vwrAq3ov5FYAojI0qz6wslzXfRbaZwTQzscLeq/bGgBqAziQ3/HlFf165a8ATpL8JsOqIvn5yO58FOHPRxmlVAn9uSuArtD6CWwD0EcvlvXzkfa56QNgq96Sc38c3cPPkQ9ovUlPQ7uO8ZGj43HA+68JrfeoH4DAtHMA7TrNFgBnAGwGUNLRsebhOVgErQkwBdr1qzeye//QeqP+oH9ejgNo7uj48+l8LNTfr7/+j6hChvIf6efjFIDujo7fzueiLbSmcn8Ax/RHj6L6+cjhfBTVz4c3gKP6+w4AMF5fXhPal5VgAMsAOOvLXfTXwfr6mvaIQ0ZsE0IIIQqpotycLoQQQhRqksSFEEKIQkqSuBBCCFFISRIXQgghCilJ4kIIIUQhJUlciCJMKVVCKfWW/ryiUmr5nbYRQhQccouZEEWYPgb2apINHR2LEOLume5cRAjxAJsK4CF9EoczAOqRbKiUeh3a7Evu0Eba+gralL2vAEgC0INkhFLqIWgDnJQBEA/gTZJB+f0mhCiqpDldiKLtfQAh1CZxeDfLuoYAngPQAsBnAOJJNgGwD8Crepk5AP5HshmAsQBm5UfQQgiN1MSFENnZRm3e6BilVBSAf/XlxwF467NZtQGwLMO0yM75H6YQRZckcSFEdpIyPLdmeG2F9r/DAG3uZJ98jksIoZPmdCGKthgAxe5lQ2pzSZ9TSvUFtFmulFKN7RmcECJnksSFKMJI3gSwRykVAODLe9jFSwDeUEqlzYTXy57xCSFyJreYCSGEEIWU1MSFEEKIQkqSuBBCCFFISRIXQgghCilJ4kIIIUQhJUlcCCGEKKQkiQshhBCFlCRxIYQQopCSJC6EEEIUUv8Pm73+TU/ZCGMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize=(8,3))\n",
    "ax = plt.subplot(111)\n",
    "line, _ = ax.plot((sol1-true_sol).abs().mean(1), c='black')\n",
    "line.set_label('torchdyn')\n",
    "line, _ = ax.plot((sol2-true_sol).abs().mean(1), c='red')\n",
    "line.set_label('torchdiffeq')\n",
    "ax.set_ylabel('error')\n",
    "ax.set_xlabel('time')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Adaptive-step bench"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.18353033065795898\n",
      "0.032141685485839844\n"
     ]
    }
   ],
   "source": [
    "t0 = time.time()\n",
    "t_eval, sol1 = odeint(f, x, t_span, solver='dopri5', interpolator='4th', atol=1e-4, rtol=1e-4)\n",
    "t_end1 = time.time() - t0\n",
    "print(t_end1)\n",
    "\n",
    "t0 = time.time()\n",
    "sol2 = torchdiffeq.odeint(f, x, t_span, method='dopri5', atol=1e-4, rtol=1e-4)\n",
    "t_end2 = time.time() - t0\n",
    "print(t_end2)\n",
    "\n",
    "true_sol = torchdiffeq.odeint(f, x, t_span, method='dopri5', atol=1e-9, rtol=1e-9)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f7806d79df0>"
      ]
     },
     "execution_count": 71,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAegAAADbCAYAAACmyQK7AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAABMAUlEQVR4nO3dd3hUVfrA8e9JD6RCQkICIXQILUAEIl2aYgNE0UVdC+K6riuuArqsqLsrKrKWn4siYsFFRZpYqNJ7Tegghl5DeiV1zu+PMxlCTyCTSXk/z3OfJFPunLm5M+897T1Ka40QQgghKhYnRxdACCGEEJeTAC2EEEJUQBKghRBCiApIArQQQghRAUmAFkIIISogCdBCCCFEBVThArRS6gul1Dml1J4y2l+hUmqHdfupLPYphBBC2JuqaPOglVI9gEzga6116zLYX6bW2uvmSyaEEEKUnwpXg9ZarwGSi9+mlGqslFqslNqulFqrlGrhoOIJIYQQ5aLCBeirmAo8p7XuCLwEfFyK53oopbYppTYppQbZpXRCCCFEGXNxdAGuRynlBdwKzFZKFd3sbr1vCPDPKzztlNZ6gPX3BlrrU0qpRsAKpdRurfUhe5dbCCGEuBkVPkBjavmpWuvIS+/QWs8D5l3ryVrrU9afh5VSq4D2gARoIYQQFVqFb+LWWqcDR5RS9wMoo11JnquU8ldKFdW2A4CuwD67FVYIIYQoIxUuQCulvgM2As2VUieVUk8Cw4EnlVI7gb3AvSXcXUtgm/V5K4G3tdYSoIUQQlR4FW6alRBCCCEqYA1aCCGEEBKghRBCiAqpQo3iDggI0OHh4Y4uhhBCCFEutm/fnqi1DrzSfXYN0EqpF4ARgAZ2A49rrXOu9vjw8HC2bdtmzyIJIYQQFYZS6tjV7rNbE7dSKhT4KxBlzantDDxor9cTQgghqhJ790G7AJ5KKRegBnDazq8nhBBCVAl2C9DWDF6TgOPAGSBNa7300scppUZac2VvS0hIsFdxhBBCiErFbn3QSil/TEKRhkAqJpf2w1rrGcUfp7WeilkMg6ioqMsmZefn53Py5Elycq7adS1ukoeHB/Xq1cPV1dXRRRFCCGFlz0FifYEjWusEAKXUPMyiFzOu+axLnDx5Em9vb8LDwym2WIYoI1prkpKSOHnyJA0bNnR0cYQQFVVyMnz1FRw8CKGh8PLLIBf1dmXPAH0c6KKUqgGcB/oApR6inZOTI8HZjpRS1K5dG+leEEJc1fTp8OyzkJUFgYGQkABbt8KsWeDh4ejSVVn27IPeDMwBYjBTrJywNmWXlgRn+5LjK4S4Iq1NTfmxx+CWW2DnTjh3Dj7+GH7+Gf7xD0eXsEqz6yhurfVrWusWWuvWWutHtNa59nw9e0hNTeXjjz8uk32Fh4eTmJhYZo8TQgi7+uc/4Z134Omn4ddfoW1bc/szz5jb3n/f1KSFXUiqz+sobYAuKCiwY2mEEKKcfPMNvP46/PGPpsbsckmP6DvvQHAwPPQQxMc7pIhVnQTo63j55Zc5dOgQkZGRjB49mtGjR9O6dWvatGnD999/D8CqVavo3r0799xzDxERERQWFvLSSy/RunVr2rZty0cffWTb30cffUSHDh1o06YNBw4cACApKYn+/fvTqlUrRowYQdEKY+PHj+eDDz6wPXfcuHF8+OGHrFq1il69ejF06FBatGjB8OHDkVXJhBBlZscOePJJ6NkTpk4FpyuECl9fmDsXzpyBO+6A9HT7liknB86ft+9rVDAVKhf39YwaNYodO3aU6T4jIyMvCoKXevvtt9mzZw87duxg7ty5TJkyhZ07d5KYmMgtt9xCjx49AIiJiWHPnj00bNiQTz75hKNHj7Jjxw5cXFxITk627S8gIICYmBg+/vhjJk2axLRp03jjjTfo1q0b48ePZ8GCBXz++ecAPPHEEwwZMoRRo0ZhsViYOXMmW7ZsYffu3cTGxrJ3715CQkLo2rUr69evp1u3bmV6bIQQ1dD58/CHP0Dt2jB7Nri5Xf2xXbrAnDlw990weDAsXAju7mVTjuPHzSC0JUtg2zZITTW3+/tD+/YwYAAMH25GlFdRUoMuhXXr1vHQQw/h7OxMUFAQPXv2ZKu1/6VTp062aUrLli3j6aefxsXaJFSrVi3bPoYMGQJAx44dOXr0KABr1qzh4YcfBuDOO+/E398fMH3RtWvXJjY2lqVLl9K+fXtq165te7169erh5OREZGSkbV9CCHFTxo2D/fvNlKrAK67hcLE77oAvvoAVK+Dee81I75uRkACjRkHTpjB6tGk+HzYMJkww2wMPmClfY8dCeDiMGGFq8VVQpapBX6um62g1a9Ys0ePcrVeXzs7OJeqvHjFiBF999RVnz57liSeeuGw/pdmXEEJc044d8OGH8Kc/Qb9+JX/eo49Cfj6MHGme98svUKxiUiIZGfDeezBpEmRnw+OPw9//Do0aXfnxcXHwf/9nmuDnzoXJk03NvwqRGvR1eHt7k5GRAUD37t35/vvvKSwsJCEhgTVr1tCpU6fLntOvXz8+/fRTW9As3sR9JT169ODbb78FYNGiRaSkpNjuGzx4MIsXL2br1q0MGDCgrN6WEEJczGIxc51r1zY11dJ68knTJL59O/ToAdYxNteVmwsffQSNG5tBaf37w549MG3a1YMzQJMmJkDv2gWtWpnm7hdegCpUWZEAfR21a9ema9eutG7dmo0bN9K2bVvatWvHbbfdxsSJEwkODr7sOSNGjCAsLMz22KLgezWvvfYaa9asoVWrVsybN4+wsDDbfW5ubvTu3ZsHHngAZ2fnMn9/QggBmGQkGzbAxImmn/dGDBkCixbB2bOmn/gf/zC/X0leHkyZYgLtX/9qguymTaY23LJlyV+zWTNYuRKefx4++MDU4K9TKaosVEUa/RsVFaUvXQ96//79tCzNP6uKsVgsdOjQgdmzZ9O0aVO7vU51P85CVGvJydC8uQl2a9deedR2aZw5YwLmnDnm786doXVrE4y9vWHvXhOI4+MhOhreeAP69oWbTZr0v/+ZPulOncy87UqQ5UwptV1rHXWl+ypVH3R1s2/fPu666y4GDx5s1+AshKjmJkwwQfrjj28+OAPUrWtGYB88CDNnmmD5008mCxmApycMHHihz7qsshk+8ogZdf7ggyb72bffls37cRAJ0BVYREQEhw8fdnQxhBBV2cmT8N//muDWrl3Z7rtZMxg/3mwAaWlmlHdwsP0C57BhcOzYhVHeb79tn9cpBxKghRCiOvvXv8wAsddft/9r+fqazd5Gj4ajR022s/BwMyq9EpIALYQQ1VVcHHz+ucmtHR7u6NKUHaXMCO/jx83I9Pr14c47HV2qUqu8jfNCCCFuzvjxJvNXVVyVysXF9H9HRppm7+3bHV2iUpMALYQQ1dHvv5sA9vzzEBTk6NLYh5eXSZpSu7bJeFbSudkVhATo63D0cpO33nqr7fbRo0fTqlUrRo8eTUJCAp07d6Z9+/asXbu2TMonhKhGpkwBZ2czB7kqq1sXli41zd59+kAlGngrAfo6HL3c5IYNG2y/T506lV27dvHuu++yfPly2rRpQ2xsLN27dy/T1xRCVHHnz8OXX5oFLq6QbKnKad4cli0zK2L16QMnTji6RCUiAfo6HLncJICXlxcA99xzD5mZmXTs2JF33nmHMWPG8OOPPxIZGcn58+dZunQp0dHRdOjQgfvvv5/MzEwAFi9eTIsWLejQoQN//etfueuuu8rr0AkhKqrvv4eUFPjznx1dkvLTpo1ZGSs52SRFqQRrWFeuUdyjRplk7mUpMtKkh7sKRy43WdxPP/2El5eXbbnNoKAgtm3bxn//+18SExP597//zbJly6hZsybvvPMO7733HmPGjOGpp55ixYoVNGnShGHDhpXlkRNCVFaffGLSafbs6eiSlK+oKLMkZv/+0KsX/PijmatdQUkNuhTKe7nJktq0aRP79u2ja9euREZGMn36dI4dO8aBAwdo2LAhTZs2RSllew0hRDUWEwNbtpi5wWWVwasy6doVFi+GxES45RaTezw729GluqLKVYOuhstNloTWmn79+vHdd99ddPuOsm5tEEJUfp98AjVqmCUiq6vu3c20q2eeMRnHvvjC5PG+5RZHl+wiUoO+DkcvN1kSXbp0Yf369cTFxQGQlZXFwYMHadGiBUePHuXQoUMAlwVwIUQ1k5oK33xj1k3283N0aRwrLAwWLDB5wjMzzQIbXbuawXPp6Zc/vrAQ/vlPkxq1nEiAvg5HLzdZEoGBgXz11Vc89NBDtG3blujoaA4cOICHhwdTp07lzjvvpEOHDtSpU6dU+xVCVDFff21GcD/zjKNLUnH07Qu7d8O770JSEjzxBNSqBbfeCq++CqtXw9atcO+98NprZoBdOZHlJquRVatWMWnSJH755ZfL7pPjLEQVpzVERJhc2Js2Obo0FZPWsHGjWdN62TLTV2+xmPucneGjj8r84kaWmxRCiOpu3TqTSevLLx1dkopLKVNzvvVWs4hIWhqsWQP5+aYJvF69ci2OBOhqpFevXvTq1cvRxRBCOMKXX5rUl/ff7+iSVB6+vnD33Q57eemDFkKIqi4zE2bNggcegBLOOBGOVykCdEXqJ6+K5PgKUcXNnQtZWfD4444uiSiFCh+gPTw8SEpKkiBiJ1prkpKS8PDwcHRRhBD28tVX0KSJmUYkKo0K3wddr149Tp48SUJCgqOLUmV5eHhQr5wHPwghysnhw7BqFfz739Uzc1glZtcArZTyA6YBrQENPKG13liafbi6utpSaAohhCilr782gbk6Zw6rpOxdg/4QWKy1HqqUcgNq2Pn1hBBCFNHaBOg+faB+fUeXRpSS3fqglVK+QA/gcwCtdZ7WOtVeryeEEOISsbFw5Ag8+KCjSyJugD0HiTUEEoAvlVKxSqlpSqnLxvcrpUYqpbYppbZJP7MQQpShuXNNBqx773V0ScQNsGeAdgE6AJ9ordsDWcDLlz5Iaz1Vax2ltY4KDAy0Y3GEEKIa0RrmzDHrHgcEOLo04gbYM0CfBE5qrTdb/56DCdhCCCHsbe9eOHgQ7rvP0SURN8huAVprfRY4oZRqbr2pD7DPXq8nhBCimLlzzejtwYMdXRJxg+w9ivs54BvrCO7DgKSxEUKI8jBnDnTrBldYEldUDnYN0FrrHcAVl9ESQghhJwcPwp498MEHji6JuAkVPtWnEEKIUpo71/wcMsSx5RA3RQK0EEJUNb/8AlFRkpykkpMALYQQVUlSEmzaBHfe6eiSiJskAVoIIaqSJUvAYoGBAx1dEnGTJEALIURVsmABBAaaJm5RqUmAFkKIqqKwEBYvhjvuACf5eq/s5D8ohBBVxZYtkJwszdtVhARoIYSoKhYsMItj9O/v6JKIMiABWgghqoqFC+HWW8Hf39ElEWVAArQQQlQFp0+b9Z9lelWVIQFaCCGqgkWLzE/pf64yJEALIURVsHAh1KsHrVs7uiSijEiAFkKIyi4vD3791TRvK+Xo0ogyIgFaCCEqu3XrICNDmrerGAnQQghR2S1cCG5u0KePo0siypAEaCGEqOyWL4du3aBmTUeXRJQhCdBCCFGZJSfDzp3Qu7ejSyLKmARoIYSozNasAa2hVy9Hl0SUMQnQQghRma1aBZ6ecMstji6JKGMSoIUQojJbtcqk93R3d3RJRBmTAC2EEJVVcjLs2iX9z1WUBGghhKisVq+W/ucqTAK0EEJUVtL/XKVdN0Aro355FEYIIUQprFoFXbuaJCWiyrlugNZaa2BhOZRFCCFESSUmSv9zFVfSJu4YpZS0oQghREWxZo35Kf3PVZZLCR/XGRiulDoGZAEKU7lua7eSCSGEuLpVq6BGDYiKcnRJhJ2UNEAPsGsphBBClM6qVSb/tvQ/V1klauLWWh8D/IC7rZuf9bbrUko5K6VilVK/3HAphRBCXJCQALt3S/N2FVeiAK2Ueh74Bqhj3WYopZ4r4Ws8D+y/seIJIYS4jPQ/VwslHST2JNBZaz1eaz0e6AI8db0nKaXqAXcC0268iEIIIS6yapVZWlL6n6u0kgZoBRQW+7vQetv1fACMASylK5YQQoirWrnS9D+7ujq6JMKOSjpI7Etgs1LqB+vfg4DPr/UEpdRdwDmt9XalVK9rPG4kMBIgLCyshMURQohq6tw52LsXHn7Y0SURdnbdAK2UcgI2AauAbtabH9dax17nqV2Be5RSAwEPwEcpNUNrfdFZpbWeCkwFiIqK0qUrvhBC2FFKCvz8M2zdapqUo6NhwADw8HBcmaT/udq4boDWWluUUpO11u2BmJLuWGv9CvAKgLUG/dKlwVkIISqk/Hx4+214913IyABvb8jJMbcHBZn7/vhHUCXp6StjRf3PHTuW/2uLclXSPujlSqn7lHLE2SiEEOUoORluvx3Gj4d+/WDzZkhLg6wsWLwYGjaExx+HJ5+EvLzyL9/KldC9u/Q/VwMlDdBPA7OBXKVUulIqQymVXtIX0Vqv0lrfdUMlFEKI8pKdDXfeCevWwVdfwdy50KmTqSm7uprm7fXr4bXX4Msv4aGHoLDwurstM/HxsG8f9OxZfq8pHKYkq1k5AbdrrZ201m5aax+ttbfW2qccyieEEOVDazPwassW+O4704R9JU5O8Prr8N57MG8ejB1bfmVcscL87NOn/F5TOExJVrOyAP8th7IIIYTjTJkCP/wAEyfCkCHXf/yoUfDss/Cf/8Av5ZQocfly8PODDh3K5/WEQ0kftBBCxMXB3/5m+p7/9reSPUcpE5zbtoUnnjDTn+xt+XIzetvZ2f6vJRyuNH3Qs7jBPmghhKiwtIbnnjOLTnzxRelGZru7wzffQGoqjB5ttyICcPgwHD0Kffva93VEhVHSAO0LPAb829r33AroZ69CCSFEuZk/34zO/uc/oW7d0j+/dWsTnL/+2kyBspfly81P6X+uNpTW188NopT6BJOu8zatdUullD+wVGt9S1kWJioqSm/btq0sdymEEFdXUGACrLMz7NwJLiVNrniJ7Gxo1Qp8fSEmxgwkK2sPPghr18LJk46Zfy3sQim1XWt9xaTqJT2LOmutnwVyALTWKYAsQiqEqNymT4fffoM337zx4AxQowZMmGCC/HfflV35ilgspgbdp48E52qkpAE6XynlDGgApVQgsgCGEKIyy883zdqdOsG99978/oYNg/bt4R//MPsuS7t3Q2KiNG9XMyUN0P8H/ADUUUq9CawDJtitVEIIYW8zZ8Lx4yZjWFnUSp2c4F//MgO5yroWLf3P1VKJ+qABlFItgD6YZSaXa633l3VhpA9aCFEuLBYzPUop2LWr7JqNtYZ27Uzf9p49ZdcXfeedcOgQHDhQNvsTFUZZ9EGjtT6gtZ6stf6vPYKzEEKUmwULzJKNY8eWbZ+uUvDyy7B/v1kFqyzk55sVrKT2XO3YYaihEEJUcO+8Aw0amH7jsvbAA2ZBjbfeMjXqm7VlC2RmSoCuhiRACyGql3XrzIIXL75onxWhXFzMvOjNm2H16pvf37JlpmYu6z9XOxKghRDVy4cfQq1aJj2nvTz+uFk3+t13b35fy5eb3Nu1at38vkSlIgFaCFF9JCbCjz/Co49CzZr2ex0PD3jqKZOh7MSJG99PVhZs2iTpPaspCdBCiOrj22/NoCt71p6LPPGEGS3+5Zc3vo+VK015JUBXSxKghRDVg9bw+ecQFQVt2tj/9Ro2NIH188+hsPDG9rF4sclS1r172ZZNVAoSoIUQ1UNsrJnz/Pjj5feaI0aYZChFiUZKQ2tYtAhuu82smiWqHQnQQojq4csvTaB76KHye81Bg6B2bZg2rfTP/f13s8TkHXeUebFE5SABWghR9eXkmHWbhwwBf//ye113dzMgbf58SEgo3XMXLTI/JUBXWxKghRBV348/QkpK+QwOu9STT5qBXp9/XrrnLVoEzZubvmxRLUmAFkJUfV9+CWFhpj+3vLVqBf37wwcfwPnzJXvO+fMmycntt9u1aKJikwAthKjaTpyApUvhscfKbvGK0nrlFYiPN+tPl8SqVaZZXpq3qzUJ0EKIqm36dDMi+rHHHFeGnj0hOhpeew3OnLn+4xctAk9P8zxRbUmAFkJUXUWJQm67zbF9uUqZkdyZmfDggyZD2NXk55s+8169TEYyUW1JgBZCVF1r15qpSuU59/lqIiJg6lRTpvbtzSpVVzJjhpk7/ec/l2/5RIUjAVoIUXV98QX4+JjpVRXB8OGwYoXpX771Vhg/HrKzL9yflwdvvgkdO8KddzqunKJCcHF0AYQQwi7S02H2bDMPuUYNR5fmgl69TEazv/wF/vUv0wT/xBMQGgqffAKHDsEvv5hmcVGtSQ1aCFE1zZplpis5Yu7z9fj5mabs1avNXOd//QuefhpOnzZJTaT2LLBjDVopVR/4GggCNDBVa/2hvV5P2FleHuzZA9u3m5zGv/0GcXGQnAx//KNZWq9tW7nqFxXHF1+Yft9bbnF0Sa6uRw9YtsxkGcvNNWtIu7o6ulSigrBnE3cB8KLWOkYp5Q1sV0r9qrXeZ8fXFGUpNdXUQtatg59/Nn8D+PqaL76ePc0o2alTYfJk02+2erV919kVoiQOHICNG2HSpMpx0RgY6OgSVG7nz5uKw+bNsG8fdOtmLn6cnc3cdycncHExqVfd3cHNzYyWT0yEpCRo1Ai8vK6839OnYf9+iIkxFZSRI8uthcNuAVprfQY4Y/09Qym1HwgFJEBXBmvXwsMPm9GktWvDPfeYkzIqykxXKf6lN2mSCeSjRsHf/gaffuqwYgsBmHPQ1dWcw6LqiY83a2WvXWuC8s6dUFBg7vPxKf3iJD4+0K+facXIyIBz58x89aJKCZjvvObNIS2tzN7G9ZTLIDGlVDjQHthcHq8nbkJmJowbBx99ZK4q1683CRauVQsJDoa//hVOnYKJE+HgQRPQQ0PNaNVBg8wH4FK5ubBpk/ki9fAwH4gzZ8wVbUqK2VJTzdq9jz9u9ifE9WRlmYFXQ4eaJmNRNezZA99/Dz/8AHv3mtu8vEwXxksvQefOZgsONt9bhw6ZdbgtFrPl55vvnKLNxcW0XPj6wk8/wdat4O1ttpYtoU8fCAmBunWhcWOIjDT3lSOltbbvCyjlBawG3tRaz7vC/SOBkQBhYWEdjx07ZtfyiGtYvdpkWzp2zIwwnTDhys0+V5OfD++/b5q8Dx26cHvTpjB2rPmQZGebAL57t2mSysy88r5cXc2qQ97eF/ZVr55pRo+Ohr59zQfG2flG362oqj77zDRDrl1rmjpF5ZSXZ4LmihUmMO/da5qqe/WCAQOgd28zn9ylck9GUkpt11pHXfE+ewZopZQr8AuwRGv93vUeHxUVpbdt22a38oir0Br+/W94/XVzpfjFFzf/xXbunGmGOn3aDCKLj79wn4+PqRVHRppFBFxczBVt3bpmCww0aQ6Lau0HD5o+8JgY2LbN/A1Qq5b5kHbubAao9e9fOfobhf1oDR06mIvBHTvkfKiMCgvN0qD/+IfJow7QvTsMG1YlW0UcEqCVUgqYDiRrrUeV5DkSoB0gIwOeecZ8IB5+2MzDLE2tuSSys+HsWTMww8PD9GnfzBfn2bPmqnrZMrMVfYiffx7ee89xCyIIx9uwAbp2NX3QI0c6ujSitJYuhTFjTJ9yVJRpeevevcoF5eKuFaDt+U3WFXgEuE0ptcO6DbTj64nSmjfPND9/843JXvT112UfnMEkiWjUyDRRBwTcfK0mOBj+8AdT0z9+3PRV//Wv8OGH5kM9fbqpkYvqZ/Jk00IzfLijSyJKY8cO0wI2YIBJMDNzphn8VQVrzKVhz1Hc6wBpX7pZv/9umnWdnc3owdRUU+vNyTFTAE6cMNMGJk2C+vUvf35WlmkabtvWPC8uzjxn61YzECwqyiTm79y53N9amfHzM2vtRkaa4/DYY+Yq/A9/gLvuMu/Rx+f6FwZ5eRemY1QGhYWmxrhhg+mnT042m1LQoAE895zpo6suTp82mcOeeUam+lUWx4/Dq6/C//5nxpy8/775/7m7O7pkFYLdB4mVhjRxF5OdbfL0fvjhhekDRZQyTcXu7maU4fHjJqgMHGiCVMuWpjl5zRozYCsh4cqv8cc/mqbAqvRh0BqWLzdN9b/8YoIumCD+6qtw//3meB07ZhZROHLE/Dx8GE6eNIPS7rjD9GO2aWN+1qnj0Ld0mdRU+O9/zVbUtx8UZLoO/P3N37t2mQuyAQPMWsQ9epRff2zRiNn8fHNulVfijRdfNJ+XgwdNi42ouFJT4a23zP8LTPfUK6+Yz2k147BBYqUlAdpq3TozrSguDkaMMM23Tk6mFujvb2oHxb9s4+LMyb1584X+WDDP6d/f7OvgQfPc5s0hLMw0N1ek/MT2kJ5uplvs2WPmTC5adPlj6tY1X+YNG5rt5ElYssTUxopERppAN2CAWeDAURc0iYmmpeCjj8x7GzjQXGT16WOCc3GpqeYi5YMPzIC9Ll3MOXLXXWXfR793L8yZY1pq9u41FztF3yu+vuYcfvBBc7Fjr/EBiYmm1WDIEFMbExWP1ub/9M03JrVpSgo88oj5PSzM0aVzmGsFaLTWFWbr2LGjrtYKC7V+5RWtldK6YUOtV64s/T4SE7XetEnrFSvM78KwWLSeP1/rTz/VetEirffu1Tor6+qPT0rSevVqrSdM0LpnT61dXLQGrWvW1Pqee7T+7bdyK7q2WLT+73/Nayul9dChWsfElOy52dlaT56sdXi4KX+PHlofOHDzZcrI0Pq997Ru3drs18lJ64gIre+/X+tx48xxmzhR6wce0NrZ2TymUSOt//1vrU+evPnXv9Rzz5ljs3dv2e9b3LxVq7Ru3NicB6B1v35ax8Y6ulQVArBNXyUmOjwoF9+qdYBOStJ60CDzLxkxwnwBioojPV3rH3/U+s9/1rpWLa19fbWeMUPrggL7vm5amglyoPXtt2u9b9+N7Scvz1yc+Plp7eam9T//qXVubun388UXJhB7eJgyRUebi4ezZ6/+nIQErb/6SuvbbrsQzO+6y1ww5eXd2PspbvNmE5yfffbm9yXK1vnzWr/4ovn/NGmi9fvva71mjaNLZZOcnKzfffddPXnyZH3ixAmHlEECdEW3aZPWdeuaWtqHHzq6NOJ6jhzRun178/Fp1swErbIINJfasUPrpk1NDfTtt00Ly806c+ZCwI+I0HrxYlNDv56MDK2feMI8r1MnrV94QeuNG0v/+nFxWv/97+Z8B62DgrQeM8bUsJKTS1aW4vbtM61NoaHmYkZUDBaL1j//bM5f0PpPf9I6M9PRpbJJSEjQ48aN076+vhqzmJP28PDQU6ZM0ZnlXE4J0BXZyZPmS6phw5I3WwrHKyjQetYsrdu1Mx+jBg1MU/L58ze/b4tF688+M7XUkBD71Dh+/tmUGbTu3l3ro0ev/tiNG03tRykTXMui1SA/35Th3nsvNIGD1rVraz16tNa//KL1unVab9hgashXOq6LFmnt42M+P5s333yZRNnYudO09oDWzZtrvWSJo0tkc/LkST1q1Chdo0YNrZTSQ4YM0bGxsXr//v26b9++tmAdFhamBwwYoJ977jn96quv6rFjx+oxY8boefPm6cQy7jqUAF1R5eZq3aWL6Vvcs8fRpRE3wmIxwSQ62nycAgNNU+v332t9+nTp95eZqfUjj5h99e2rdXx82Ze5SG6u1lOmaO3trbWXl9YDB2r9xhtaz5ljvlRXrTJN+s7OWoeFmT55e0hMNMfwP/8x/evFA3bR1rKluVA4edL0XT7/vGkqb9dO62PH7FMuUTp79pgxCGAunP7znxvrRiljhYWFevny5frBBx/Ubm5u2tnZWT/66KN63yXdRQUFBfrHH3/U//rXv/Tw4cN1hw4dtLe3t1ZKaVdXV+3m5qYBPXHixDIt37UCtIzidqQXXjCjbGfNMtN/ROWlNaxaBR9/bKZ35eSY0d5PPAGdOkGTJtCu3bWT7W/cCE8+aZZKfP11s2hJeeQaP3zY5F3fuNEsq1f8O8HV1byHd94xI7LLQ0KCmdedlmbKkpBgplAVny7o4gKPPmqm6dgjuU5x+flmVH9YmKQOBXMsfv3VTAXt0MHMHPnkE5Pdz8vrwqp2RVP+ysmBAweYO3cuK1eupHPnzgQFBbF582bWrl3LiRMn8Pf355FHHuGFF14gPDy8xPstipH5+fls3bqVBg0aUK9evTIrt0yzqojmzjVZcooyYImqIzfXTO36+GOTna1oHntoqJlCV/TlUFAA8+fDa6+ZudmZmWb621dfmalTjpCWZuaGZ2WZ8rRpY+baO1p8vEnrmpVl5spGR185MU9Z0NpcJC1bZgLRqlVmTnlgoMn/XlBg5nrXrGkuuPz9TXa74GAzba/o9+Bg87+uLIlvrmfPHpOp7csvL8/UFxZmUqv+6U+XT/mzoxMnTvDZZ58xZ84c9u/fD0BERAQHDhzAYrEQEhJCly5duO+++xgyZAgeHh7lVraSkgBd0cTFmVWZWrY0yUTc3BxdImEv+fkmKcquXaZ27ONjkoYkJ8OWLWZeaESEmWMdGGiyf9m7Rigul5trEtzMmweLF5ulU8EsHtO3L7RuDdu3m+x9Li6mJp2VZQJ3UpK5gIiPN9ndimva1Fyo9e5dOVdeS001yzt+8YW5uHR3N3OX//pXc6G0fr1pWenfv1zfX2xsLO+++y6zZs1Ca02PHj247777GDx4MKGhoZw7d468vLwyrenaiwToiqSw0KwU9dtvJv9sNZ6gX+1s3Gi+2BITzZdbRIRpRbn77qpTy6pMtDZpUr/4wiRaSU83F1ADBkC/fiYwN2xY8v0VFppgffas2Y4dM9myjhwx+42ONgt5dOtmuj0qajrS1FST1GfmTHOxkpdnLjSeftqk0S3HGvKlkpOTeeWVV/jss8/w9vZmxIgRPPfcc6Vqsq5oJEBXJBMnmhVavvvOZFcSQpSfs2fhpZdMf6mTk6kpe3mZC6X77zddC2WZKS4ry3RjrF9vaqB79pgLA2dnk9UvIsK0pEVEXFhm1d8fWrSwX3+31qZLZcUKszVpYoLuli1m++0387jQULPE40MPmRY/B/e/z5o1i2effZaUlBSef/55xo8fj295jYuwIwnQFcW+fWbxgrvvNkn9ZcCJEOWjsBCmTYOXXzaDm+67z9x2++0mMJdXt0JqqmlJKUpBu2+fGRBnsVz8uPvvNxfzDRqU7HsiN9cE3kv7WPPzzSCuw4dNGtb1602rQVETfkCAqfVrbfrMO3c2tfsePSjo1InYnTvZuXMnJ06cIDs7m+zsbNzc3GjcuDGPPvooPj4+ZXJYriUvL4+XXnqJjz76iE6dOjF16lTatWtn99ctLxKgK4KCAtPEdfSo+aBUtAUYhKiKtIYFC+Dvf4fdu6FXL5gyxdReK4qcHLNqXUqK6ePevNnkpy4oMP27rVubrVkz0yUWGmruS0szNeHFi03/uZOTyc2el3dhAZjjxy/uF2/QwOSTL2pqb9vWBOicHLNfpdixYwdTpkxh5syZpKWl2Z5ao0YNPD09yc3NJTMzEz8/Px5++GFGjhxJmzZt7HJoYmNj+ctf/sKGDRsYNWoUEydOxLW8Fl8pJ5KLuyJ46SUzP3D2bEeXRIjqYc0arbt2NZ+7Jk20njmz9JnKHGXfPq0//tjMQ+/eXWt//8vnhhdt4eFm7v2wYWZueJ06Jr/CH/5g8qJ//rnJ63/q1DVfct26dbpr1662rFqPPPKInjlzpj506JAuvCSL3ZYtW/SwYcO0u7u7BvR9992nd+/eXWZvPz4+Xt99990a0N7e3vr7778vs31XNEiiEgebN88c6j//2dElEaJqs1i0XrhQ6969zWeubl2TjMUeqVjLk8ViUqHu3Kn1ggVa//qryZ525MjFFx2lzPKWnp6uN27cqB944AEN6JCQEP3+++/r5OTkEj0/KSlJv/rqq7aEHk888YROSEgoVRkutWzZMh0cHKzd3d31hAkTdEpKyk3tr6KTAO1IWVla16tncjdXgKw6QlQ6V6v1njtnFizZvt1kF5swweQXB5Obe9Kka69YVs3NmzdP16xZUwPa09NTv/baazechzoxMVG/+OKL2sXFRdeuXVtPmzZNW0rZWpGXl6f//ve/a6WUbtGihd6xY8cNlaWykQDtSG++aQ6zvdIkClEVHD9uFh154QWzMMiYMWZ1t4gIrT09TRrS1au1/u47rV97zayMVbQEaPEtOlrr6dPlYvgaLBaLfuedd7RSSnfu3FnPnz9fn73WamSlsHv3bt2tWzcN6N69e1+WTvNqjhw5oqOjozWgn3zyyXJfsMKRrhWgZZCYPZ07Z6Yw9OljJvsLUZ1kZZlBTwEBF25LT78wgGnfPpPAJTbWJO8BM8UpN9ck72nSxMy/DQ6GGTPM/sCMam7d2ozAHjLEDLBSygwAqwSJKcqaxWJh69ateHh4XHd0c3Z2NiNGjOC7777jgQce4Msvv6RGjRplWh6tNZ9//jkvvvgimZmZ3HnnnbRr1w4vLy969+5Np06dLirPZ599xmuvvYbWmk8//ZQHq9n0UxnF7SjPPgtTp5pR282aObo0oirS2mQlO37cJMY4dsxkvEpPNyN6GzUyc26dnc0o36Lfi6bfHDtmAuhdd5n7ryU//8IUnSFDLp7So7VJwHLwoEnAs2CBmWObnw/du5tAffiweUxxjRubkcTdu1/I2JWZCTVqXJyZ6sQJ874aNzaB29OzzA5hZaS1JjY2lhkzZjBz5kzOnDmDl5cX+/fvv2r2rAMHDvDQQw+xc+dOJkyYwNixY1F2nOqZkJDAO++8w08//URcXBxFsaZ169Y4OTmRk5PD4cOHKSgooF+/fkyZMoVGjRrZrTwVlQRoR9i/3+QxfuYZ+OgjR5dGVFbnzsHSpbBwobnQ69zZBNJjxy4E5aKaZZGgIBPcTp8u+es0aWKm2RQP4kWBPCvLTOmJi7uQg7l+fZNnOiXFzO1NSbm4HI0bm/n+NWqYsgcEmIuFxo3Nz0aNTO34WouHiIsUFhYSExPDvHnzmDt3Lr///juurq4MHDiQ22+/nb/97W/079+fDz/8EG9vb1xdXUlNTeXAgQP89NNPTJs2jRo1avC///2PgQMHlnvZMzIy+PTTT1mzZg1ubm64urrSuHFj+vfvT69evex6sVCRSYAub1qb3Lu7dpmsPIGBji6RKE/x8WYxlORkM2/1+HET3Jo0gVdeuVAzzMkxc1i/+MIssvDQQyYt5J49Zs7unj1w5ox5bJ06pqa5ZYtp/g0LM3NaGzS48HtYmNkCA8082aJFHgoLzWaxXPjd2dk0BzdoYObdfv21qeUW3V+0ubmZIOvlZVqBOnc2f3/wgZlv6+dnMl/5+Zl9NWtmsmCFh0sinhuQl5fH3r17CQwMJCUlhe3btxMTE0NMTAw7duwgKysLFxcXevfuzdChQxk6dCi1atUC4K233uLvf//7Fffr4eHBfffdx6RJkwgODi7PtySuQwJ0efv6a9O8OHUqPPWUo0sjylpRgAOTlWrfPhNQN2826Rz37r38OUFBJnDfdZcJxgcOQEyMqXWGhJhkEUW1Uw8PaNXKNPe2amUu9jp0uH4TtLiu7du3M378eGJjY2nevDkJCQkkJSVhsVgICgoiODiYZs2a0bVrV7p160Z9e62YZXXu3Dm2bNnC1q1b2bp1K+vXryc9Pf2ix9SsWZP27dvTsWNHbrnlFu644w5bUC7OYrGwePFizpw5Q0ZGBgUFBfj4+NCkSROioqLKJeuXKD0J0OUpKcnUIJo2NV/W8qVa+eXnm9aQjRvhp59M1qamTU2NNC7uwvrJPj4mS1O3bnDvvaYWeeqUaQ6uUcOkbhw3ztQ2W7Qw60PfeadZmOHcOROwW7QwCzRUxpWPKrA9e/Ywfvx4fvjhB2rVqsWAAQM4cuQIQUFB1LFm9Tt37hxnzpxh7969ZFmb65s3b87AgQMZNGgQPXr0uOlyJCYmsmTJEpYsWcL69es5fPgwAE5OTkRERBAdHU3v3r1JTU3Fx8eHjh070rRpU5zlfKiyJECXF61hxAiYPt182bZt6+gSiZLIzDTNwTNnmlps166maTkuzmy//26af8EEz3vuMQOeXF3NOIOirXHj6wfWvDxZXrQcnT59mpdffpkZM2bg7e3Niy++yKhRo65ZmywoKGDnzp2sXbuWxYsXs3LlSvLy8njmmWd4//33cS/hYhpF363p6enMnj2br7/+mnXr1qG1JjAwkB49etC5c2c6d+5Mhw4d8KoEy4wW9YNnZWXRs2fPattvXJYkQJeHtDR4/HEznWrMGHjnnZI/d+lSMzp2xAizqsyBA6aZdNcuU1srLIRBg0w/pLu72VxcTG19xw5z/5tvmqbT8nT+vCm3v78pS0CAWQykiNamT7VBgysHroIC2LbN1Br9/C6/X2sTIAMCzGtcKjfX9Ml++61pYn7zTejZ09yXlWVGFB88aGq+W7aYkcJhYSb4Fi0JePiweUxhIdSqZZqXT582QbRRowtTfTp3NrnU69eXvtVKICcnh/fee48JEyaQn5/PqFGjGDt27BWbhq8nKyuLN954g3fffZegoCBGjhzJ008/TWho6BUfn5qayrRp0/joo484ffo0Siny8/Np3rw5Dz30EAMHDqRjx444VfDWtezsbPbs2cOOHTts265du2ytCwMHDqRXr164u7vj7u6Oh4cHHh4eF/3u4eFBbm4uZ86cIT4+nn79+hEREeHgd1axSIAuD48/buZqTpgAf/vb9WtS2qxFW/jWWzgvWACAxcUFJ4vlwso2rq7QpYsJHhs2XHk//v4mUHp5mfnWDRuawBIeboLM8eMmWD35pNnfzbJYTF/rDz/AV19BQsKF+5yczIj1li1hzRqzpOZvv5km3PHjTX9terop044dsGyZCZYhIWad5Px8c396umkajo01P2vUgDvuMK+dmWkGPiUlmeBfUGCCau3aZl/+/ub+vLwL5fLwMDXc2FjzeCcn0yccHGwuHlq3hh49THB3cTF9xXXqVPlm5oMHDzJ//nzS0tLo2bMnHh4eFBYW4uLiQmBgIHXq1MHf379S1ZK01syfP58XX3yRI0eOMHjwYCZNmlQm03eWL1/O+++/z8KFCwFo27Yt7dq1IzQ0FFdXV1xcXNi7dy8LFiwgMzOT3r1706VLFywWC/fddx9RUVEV6lgeOnSIb775hujoaKKjo9mxYwfbtm1j27ZtxMbGcuDAASzW7yIfHx8iIyNp3749Xbp04fTp04wfP94WrEvK09OTt99+m/bt2xMaGkpoaGiJWyRKIjs7m5iYGDp27IjnVabinT9/niNHjrBjxw7WrVvHxo0bqVOnDsOGDePxxx+/7H+ktWbjxo2sW7eOgIAAunbtSvMyXGxFAnQJaK358MMP+eijj/j4448ZMGAAYP6ZsbGxzJ49m4yMDN566y0CLx2VvX49dOvGb4MG0XTu3MuujDcsWsTpbdtoPXQo6atWoWbPpsHWrdTJzCQJmATMA54GOvbsSc8//9kEjaZNOXrqFE5OToT5+5sAl5trtoIC8PfnSFYW2TExtPrf/8xgpWPHzH2XyLnrLjzmzrU1rxYUFLBkyRISExMZMmQIXl5epKenk5KSQkJCAr/99hsrVqzgxJEjPPKHP9AuORn/1aups3EjHqmpFDo5EVOnDl+5uTFu9GhCIiLg3XfNqGQwtczoaDLat8dr2jRU0QCoIkFB0K0byZ074zd9Ok5FA6s8PU1fbmCgOQY9e8KmTbB2rQnU3t7mYsTf3zQpd+5MjLc3G7dsYWRKCq4ZGeb5fn6m9tu8uakBe3iY6UC5uaZGfoPBNyMjg5iYGKKjo3Ero6bq1NRUatSoUWb7uxqtNTt27OCHH35g3rx57LUec2dnZwqLr3hUjJubG3Xr1iU4OBgPDw/S09NtU3UaNmx4xedkZWWxZMkSpk+fztmzZxk6dCju7u6kpKSQnZ2Nq6srrq6u1KhRA19fX3x9fXF1dSUzM5O8vDyGDh2KX7EWFa01aWlpnDhxgszMTDp37nzF2ufWrVsZO3YsK1eupHXr1nzwwQf06dPnsmOwdetW8vPz6dq1a4mOm8ViIScnx5bQ4/Dhw8yYMYO1a9eyf/9+zp49azt+ISEh3HHHHfzlL38hMjKyRPu/WWfPnmXBggX4+voydOjQaz42NzeXn3/+ma+++opFixbZAnBxoaGhdOjQgfbt2xMZGUlkZCTh4eGXBa7CwkLOnz9Pbm4uubm55OTk2H4W31xcXAgNDcXNzY3HHnuM9evXX7SfwMBAQkNDqVevHn5+fmitSUxMJC4ujvr16/P222/TsGFD23mjlCItLY2UlBRSUlKIj48nJiaGLVu2sGHDBrKzs2nRogVPPPEESUlJnDt3jnPnzhEfH8+ZM2c4VbTUJuDt7U3nzp05ceIEv/32G48//jidO3e2Pf7UqVNs27aNkydP2p7zf//3fzz33HM38q+6IgnQl9i+fTsRERF4enqitWbTpk28+uqrLF++HF9fX7Kzs2nXrh2nTp3i7NmzaK1tX55BQUH07dsXf39/vLy8sKSlMXLaNHRWFi2BfvfeyyOPPEJWVhZnTp8m8+uveXb/foKBDMAbyAfWuruzs2lTUvr3p9uAAXTs2JFRo0YxY8YMnn76aSwWC7GxsWzbtg1XV1fGjh1L48aNycvLIyMjg+PHj7N582Y2b94MwLBhw+jUqROFeXk4nTlD1t697IqJ4UBKCv2B94BcV1cOBgayxdOT7+Pj2ZSZSQbg6upKYWHhRR/WlsAkFxcGFgv2mcBCYL71p29YGKmpqbRq1Yq1a9fiXFiInj+fZIuFZcnJfDJrFqtXr+a5gQN58Z57SCssJCE3lxPAzuPHWb58Obt376Z1RASvPf88eHtT6OSExWIhISGBLVu2sGjRItq0acPEiRPx8PAgOzub+Ph4jh07xqFDh1i5ciW7d+8GoEOHDgwfPtyWJs9isZCUlMSBAwfYv38/Xbp0ISQkhMTERHJzc3nllVdszW0Wi4UlS5ZQs2bNiwYDJSQksG/fPvbs2cPq1atZuHAhWVlZtGzZkqFDh1JYWEhBQYFtK/rSiomJoUGDBkyePPmyptDs7GwOHjzI5s2b+emnn1i8eDHh4eG8/fbb+Pr6XrS//Px8CgoKOH36NAkJCfTp0wcvLy/i4+OJj48nKSmJwYMH06pVK9v+s7KySEpKIiwsjMLCQjZu3GgLykePHsXJyYnu3bszZMgQBg0ahJ+fH9u3b0drjbOzM3l5eSQkJHDu3DnOnj3L6dOnOXPmDLm5ufj4+LBhwwZ8fHx4/vnnCQwMxN/fn8LCQn7//XdWrlzJsmXLyM3NJTg4mLp16xIbG2srm7u7u+04XY2Pjw+NGjXCYrGQm5vLqVOnyMzMtN0/cOBARowYwZo1a1i3bh09evTgyJEj/PDDDwQEBPDaa6/xpz/9CRcXF9vxWL9+PYsWLWLhwoUcPHgQgBdffJEuXbqQlpZ2zS0uLo7MzEzGjh3L/fffb2vCLd606+LiQlpaGtnZ2Zw8eZItW7awbds2WrduzZgxY2xluZKCggKOHz/Otm3b+PXXX9m2bRsjR460fY9kZGSQkZFBZmYmqampnD59mlOnTrF//362bt3KsWPHADPQbOPGjRdl6iosLCQuLo7Vq1fz66+/smzZMlJTUwkNDeWxxx7jqaeeYtmyZZw6dYoOHTrQsWNH6tqxq8xisXDgwAFOnTrFyZMnbT+LtvT0dJRS1KpVi/DwcFavXk1C8Va6q3BxcaFt27ZER0fTvn173njjDU6cOIGbmxt16tSxbUFBQTRq1IgmTZoQERFBmzZtcHZ2RmvNP/7xDyZMmGDbp7+/P0FBQbRr147+/ftz7733kpGRgY+Pzw11lVyNwwK0Uup24EPAGZimtX77Wo8vjwD97rvv8u2YMVhatWLosGF899137N+/n9q1a/PKiy9yS2oqU0+eJDEpidDQUOrXr09kZCQ9XFxImzOHcTExrElOJiUlhfPZ2cxSikFaM7dHDwLS0vjjrl1kac1QTI04CjgXEkLBH/9Ixpo1FHbpQp2RIwlo1gytNYcOHWLmzJns3r2b3r17s2DBAhYvXoyfnx8RERHccccdxMTEMHv27Iveh6enJx06dOD2228nPz+ft956i/z8fACUUtSrV4+uXbsydOhQatasydfDh9MjLY3eStG8WNA9X78+33TqxNkGDQjTmmZnztBk924C9u5Fe3uT8/DDJOXlkdy4MedvvRXnmjVJS0sjPj6eFStW8Ouvv3L8+HFbq0J2drat2at+/fpER0czZ86cy67UPTw86NSpE9HR0cyYMeOiq9oigYGB9OnTh8WLF5OamnrZ/Z6ennTp0oWWLVuSkpLCwoULL1q/FszFR+3atXF3dychIYHc3Fw8PDzIy8vD19eXiRMnsmvXLn788UeOHDkCwO233052djb79u0jsVjmq9DQUPr160fdunWZNWsWhw4dwsXFBWdnZ1xcXGxbUQKG2NhYCgoKCAgIwMvLC1dXVzIyMjhx4oRtAFFwcDD169fn6NGj1/0iUkpxpc+rh4cHs2fPJjs7m2XLlvHtt9+SlZVF165diYuLIz4+Hjc3N/r168fgwYO5++67sVgstgsgd3d3xo8ff9mXjtaagoICUlNTOXjwIBs2bGDu3LkcOHAAZ2dnkpOTLytLeHg43bp1o2HDhpw6dYrCwkKefvppWw0sMzOTlJQUkpOTbRda+/fv59ChQ7bR0wEBATg5OdlqZrm5uRe976Lj4OzsjLu7Ozk5OdSsWZPHHnuMAQMGkJyczO+//87u3bvZvXs3hw8fRmuNu7s7vXr1YujQoWzdupWpU6dedny9vb1ttfqirX79+mRmZvLNN99c8/9zqeDgYM6ePUtERARBQUEUFhbaLoQLCgpIT08nOTmZ5ORk2wWLt7c3DRs2ZNeuXdfct1KKhg0bEhUVRVRUFF27dmXYsGF4enrSv39/4uPjOXLkCHv37iUnJwe4cP4++OCD9O3bt1KMDE9OTmb+/PmcP3+e/Px88vPzsVgs+Pr64u/vj7+/PwEBAbRs2fKiJu38/HyysrLw9fUtVbfCsWPHcHZ2pk6dOnZv0SrikACtlHIGDgL9gJPAVuAhrfW+qz3HngE6MTGR1159lYFTpnAnsAdoA3Tv3p277riD9lu3cuv8+dTUmkOBgTRYvpzsRYs4/+uv1Ni0Ca/MTIr+zdn33ce5N9/E94UX8F+0iMMhITQ4fRpn4LyrK+4WC06FhaSFhLCuSxf+eeQI51JSGD58OAUFBRw9epSjR48SFxdHUlIScOHDPHjwYObMmYOTkxMJCQnMmzcPNzc3evfuDZgmx5o1a+Lj48Phw4eZPHkyR48e5b333qNWrVq4uLjg5uZGQUEBS5cu5a233mLr1q1Mnz6d4cOHA5B36hSnZ83i8JIltNmyhcCUlIuOVXJQEKsbNGD8sWPsiY+nY8eOODs7c+LECVuLAoCXlxdaawICAujTpw9ubm64ubnh5+dHQUEBy5cvZ/Pmzdx9990MHTqU4OBg/Pz8yM3N5dChQ8ybN49ffvmFPn368Prrr+Pt7W3LOJSQkMCePXv46aefCAkJISoqyva+ir5YDh8+zPLlyzl69ChgWjdatmxpa9IqHly9vLxstTAfHx/y8/PJycmxfdH7+PiQkpKCu7s7NWvWtDW5BQcHU7t2bQoLCzl48CArV64kOzubrl27MmrUKAoLC8nLyyM3N5ezZ88SFxfHvn372Lp1K56enrRv3x6lFNnZ2eTn56OUwtXV1VYrTkxMtAWckJAQLBYL6enpZGdnX/Q/qVWrFi1btmTbtm04OzvbmgVzc3PZvn37Fc95T09PQkJCqFu3Ln5+fpw/f57Tp09z/Phx20WUl5cXOTk5eHh44OPjQ25uLnl5ebb3dKm2bdty/vx5jh07RvPmzcnOziY3N9f2xZmSkkKB9QLQx8eHwsJCsrOzr3hhUcTHx4eWLVsSGRlJgwYN+PHHHyksLCQ4OJigoCDb/OS6dety7NgxxowZg8ViwcPDg+bNm7Nz587L9unk5ETTpk1p27Ytbdq0ISoqip49e9qaqbXW7N27F621LRB7e3tfc+BWTEwMx44du6xJt+iYeXl54ePjQ506dejYsSPBwcF8++23TJkyBYvFgrOz80Wbr68vtWrVonbt2oSHhxMZGUm7du1wcXFhzpw5HDlyBG9vb7y8vGw/fX19CQkJISgo6LJa+dKlS7n//vtxdnYmKCiI+vXr06ZNG9q0aUOXLl1o3rx5heoDF4ajAnQ08LrWeoD171cAtNZvXe05ZRWgC3JyWDtwIA3WrcM/P5/fXF2x5OcTCdQA8pXCVWvOuLjgYrHgb7HgAliADE9PfIum1FgVAnFeXiytVYuHjh+nNpAN1LQ+xwIsDw8npmNHGs6dy2ngGyDG+vwOHTrg5+fHihUrcHV1JSwsjIYNGxIeHk7Hjh3p378/DRs25L333uOll16iTZs2pKWlcfz48eLHhvDwcPLy8khPT+fw4cMcP34cZ2dnWy2xXbt2ZGZmkpiYyKFDh8jNzSU0NJSQkBC2b99OeHg4qamppKam2mq0NZ2dubewkEAgHtgEHAV8fX3p2bMn7du3Z/ny5Xh6elK/fn3q1atH/fr1ad26NZGRkcybN4/hw4fj7+9PQUEB58+ft31BN2nShC5dujBjxgxq1qxJfn4+ecUGbwUFBTFw4EC+/vrrqzZ5dunShd27d19xMIq/v7+tqTYiIoKxY8eSnZ1NcHCw7cu9SZMmREZG0rp1a9LS0nBycsLHx4d9+/bRrVs3UlJScHNzo02bNgwcOJAFCxYQExNzhZJAixYt6NWrF82aNePVV1+9YplCQ0Np0qQJffv2ZePGjbYBRcUV1ZqbNWvGbbfdxqBBg5g/fz7z58/Hz8/P9qUdGBhIs2bNaN68OSEhIbZAXvxL1mKx8PTTTzNnzhyaNGlCdHQ0ffr0ITU1lXHjxpGfn0/NmjXx8vLCy8uLunXrEhYWRnh4OB06dKBTp04cOHCAyZMnY7FYcHNzw93d/aKfPj4+NG7c2BZ0EhMTGT16NElJSZc9vlatWrRq1YqWLVvSunVrzp07x7Rp03B2dsbPzw8/Pz/8/f1tvwcEBBAUFFSqwLF//35SU1Np37497u7uzJs3j99//9026Cg0NJSwsLCrDhKqyi49P0TF56gAPRS4XWs9wvr3I0BnrfVfLnncSGAkQFhYWMei/pSbkX/+PKpGDZyBHMAT0ECOkxPnunQhbNkyUkNC8LUOdsr29qawb1/8P/gAgoOJu+02Cs+dIzU8HLp1o8lTTxEQFEReXh4/TJnCoJdfxjUnh8y6dUkfMICA0aPxaNkSgDVr1pCQkEDNmjXx9fWlefPmtqbDtLQ0vLy8rtq0pLXmrbfe4tdffyUkJIQ2bdpw++23s3PnTiZNmkRhYSFubm54eXnZmo8HDRpEQkICzz//PLm5uXh5eeHn50fTpk3p06cPvXv3pqCggHHjxpGQkGC7am/evLlt5O7ixYvx9PQkICCA2rVrExAQQHBwcIk+6FprJk+ezL59+3B3d6dGjRo0atSIzp0706pVK5RSfPPNN7bVdjw9PWnQoAGRkZG0bdsWJycn1q5dy5IlS2w1cH9/fxo0aEBUVBQBAQGcOXOGDRs22GoQPj4+1K5dmzp16tzUl1HRAJMmTZrYaiO5ubmsWbMGJycnatSogaenJ15eXtSrVw+PYotDJCQkcPr0aVuZ3dzcqF279kUrA2mtOXv2LEopW3+lm5vbNfsjhRDVS4UO0MWVZRP3pnHjiHjmGXzq1SM/IQFXf38zhaasaC3zYYUQQtyUawVoe17KnwKKJ7KtZ72tXHR5803b7672WKxCgrMQQgg7smcqm61AU6VUQ6WUG/Ag8JMdX08IIYSoMuxWg9ZaFyil/gIswUyz+kJrfYVlfoQQQghxKbuOVtFaL8TktBBCCCFEKVTsbO1CCCFENSUBWgghhKiAKlQubqVUAnDzE6EvCAASr/uo6kOOx8XkeFwgx+JicjwuJsfjYmV5PBpora841ahCBeiyppTadrX5ZdWRHI+LyfG4QI7FxeR4XEyOx8XK63hIE7cQQghRAUmAFkIIISqgqh6gp17/IdWKHI+LyfG4QI7FxeR4XEyOx8XK5XhU6T5oIYQQorKq6jVoIYQQolKqkgFaKXW7Uuo3pVScUuplR5fHEZRSR5VSu5VSO5RS26y31VJK/aqU+t3609/R5bQXpdQXSqlzSqk9xW674vtXxv9Zz5ddSqkOjiu5fVzleLyulDplPUd2KKUGFrvvFevx+E0pNcAxpbYfpVR9pdRKpdQ+pdRepdTz1tur3TlyjWNRLc8PpZSHUmqLUmqn9Xi8Yb29oVJqs/V9f29dYwKllLv17zjr/eFlVhitdZXaMHm/DwGNADdgJxDh6HI54DgcBQIuuW0i8LL195eBdxxdTju+/x5AB2DP9d4/MBBYBCigC7DZ0eUvp+PxOvDSFR4bYf3cuAMNrZ8nZ0e/hzI+HnWBDtbfvYGD1vdd7c6RaxyLanl+WP/HXtbfXYHN1v/5LOBB6+1TgGesv/8ZmGL9/UHg+7IqS1WsQXcC4rTWh7XWecBM4F4Hl6miuBeYbv19OjDIcUWxL631GiD5kpuv9v7vBb7WxibATylVt1wKWk6ucjyu5l5gptY6V2t9BIjDfK6qDK31Ga11jPX3DGA/EEo1PEeucSyupkqfH9b/cab1T1frpoHbgDnW2y89N4rOmTlAH6XKZj3iqhigQ4ETxf4+ybVPtqpKA0uVUtuVUiOttwVprc9Yfz8LBDmmaA5ztfdfnc+Zv1ibbL8o1uVRrY6HtUmyPaamVK3PkUuOBVTT80Mp5ayU2gGcA37FtBKkaq0LrA8p/p5tx8N6fxpQuyzKURUDtDC6aa07AHcAzyqlehS/U5v2mGo7hL+6v3+rT4DGQCRwBviPQ0vjAEopL2AuMEprnV78vup2jlzhWFTb80NrXai1jgTqYVoHWjiiHFUxQJ8C6hf7u571tmpFa33K+vMc8APmJIsvapaz/jznuBI6xNXef7U8Z7TW8dYvIgvwGReaKavF8VBKuWIC0jda63nWm6vlOXKlY1Hdzw8ArXUqsBKIxnRrFC3RXPw9246H9X5fIKksXr8qBuitQFPriDs3TKf9Tw4uU7lSStVUSnkX/Q70B/ZgjsMfrQ/7I/CjY0roMFd7/z8Bj1pH6nYB0oo1c1ZZl/ShDsacI2COx4PW0akNgabAlvIunz1Z+wg/B/Zrrd8rdle1O0eudiyq6/mhlApUSvlZf/cE+mH65VcCQ60Pu/TcKDpnhgIrrK0vN8/RI+bssWFGXB7E9BuMc3R5HPD+G2FGWe4E9hYdA0y/yHLgd2AZUMvRZbXjMfgO0yyXj+kvevJq7x8zanOy9XzZDUQ5uvzldDz+Z32/u6xfMnWLPX6c9Xj8Btzh6PLb4Xh0wzRf7wJ2WLeB1fEcucaxqJbnB9AWiLW+7z3AeOvtjTAXInHAbMDderuH9e846/2NyqoskklMCCGEqICqYhO3EEIIUelJgBZCCCEqIAnQQgghRAUkAVoIIYSogCRACyGEEBWQBGghqiillJ9S6s/W30OUUnOu9xwhRMUh06yEqKKseZV/0Vq3dnRZhBCl53L9hwghKqm3gcbWpP+/Ay211q2VUo9hVuKpickCNQmzNOsjQC4wUGudrJRqjEnOEQhkA09prQ+U95sQorqSJm4hqq6XgUPaJP0ffcl9rYEhwC3Am0C21ro9sBF41PqYqcBzWuuOwEvAx+VRaCGEITVoIaqnldqs/ZuhlEoDfrbevhtoa13Z6FZgdrGlbd3Lv5hCVF8SoIWonnKL/W4p9rcF873ghFn/NrKcyyWEsJImbiGqrgzA+0aeqM16wEeUUveDWfFIKdWuLAsnhLg2CdBCVFFa6yRgvVJqD/DuDexiOPCkUqpoVbR7y7J8Qohrk2lWQgghRAUkNWghhBCiApIALYQQQlRAEqCFEEKICkgCtBBCCFEBSYAWQgghKiAJ0EIIIUQFJAFaCCGEqIAkQAshhBAV0P8DBkOqdWm/VdIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize=(8,3))\n",
    "ax = plt.subplot(111)\n",
    "line, _ = ax.plot((sol1-true_sol).detach().abs().mean(1), c='black')\n",
    "line.set_label('torchdyn')\n",
    "line, _ = ax.plot((sol2-true_sol).detach().abs().mean(1), c='red')\n",
    "line.set_label('torchdiffeq')\n",
    "ax.set_ylabel('error')\n",
    "ax.set_xlabel('time')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using `ODEProblem`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 126,
   "metadata": {},
   "outputs": [],
   "source": [
    "from torchdyn.core import ODEProblem\n",
    "\n",
    "f = VanDerPol(0.5)\n",
    "t_span = torch.linspace(0, 2, 300)\n",
    "x = torch.randn(1024, 2, requires_grad=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Backsolve Adjoint"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* We are slower when the ODE is easy to solve due to the checkpointing strategy used that slows down `dt`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 127,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Your vector field does not have `nn.Parameters` to optimize.\n",
      "0.17838144302368164\n",
      "0.018263816833496094\n"
     ]
    }
   ],
   "source": [
    "prob = ODEProblem(f, sensitivity='adjoint', solver='tsit5', atol=1e-3, rtol=1e-3, atol_adjoint=1e-3, rtol_adjoint=1e-3)\n",
    "t0 = time.time()\n",
    "t_eval, sol_torchdyn = prob.odeint(x, t_span)\n",
    "t_end1 = time.time() - t0\n",
    "print(t_end1)\n",
    "\n",
    "t0 = time.time()\n",
    "sol_torchdiffeq = torchdiffeq.odeint_adjoint(f, x, t_span, method='dopri5', atol=1e-3, rtol=1e-3)\n",
    "t_end2 = time.time() - t0\n",
    "print(t_end2)\n",
    "\n",
    "true_sol = torchdiffeq.odeint_adjoint(f, x, t_span, method='dopri5', atol=1e-9, rtol=1e-9)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 128,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.6487090587615967\n",
      "0.8733053207397461\n"
     ]
    }
   ],
   "source": [
    "t0 = time.time()\n",
    "grad1 = torch.autograd.grad(sol_torchdyn[-1].sum(), x)[0]\n",
    "t_end1 = time.time() - t0\n",
    "print(t_end1)\n",
    "\n",
    "t0 = time.time()\n",
    "grad2 = torch.autograd.grad(sol_torchdiffeq[-1].sum(), x)[0]\n",
    "t_end2 = time.time() - t0\n",
    "print(t_end2)\n",
    "\n",
    "grad_true = torch.autograd.grad(true_sol[-1].sum(), x)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 131,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f78063d2130>"
      ]
     },
     "execution_count": 131,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfcAAADQCAYAAAAaqygdAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAxR0lEQVR4nO3de5wcZZ3v8c8vM5MrIUC4HQiQKFEWCIYQIiwXURbBgwur4gq65+gqi57V1y56hINHX6AcXueIh+Nlva1ZQKIICQRWAoRbiGwkQEiAkBsJJCFAuOR+TyZz6d/546marunpnq6e6Z6e6f6+X695TVd1XZ56+qn61fPUU1Xm7oiIiEjtGFTtBIiIiEh5KbiLiIjUGAV3ERGRGqPgLiIiUmMU3EVERGqMgruIiEiNaax2Asrl0EMP9bFjx1Y7GSIiIn3mhRde2Ozuh+WOr5ngPnbsWBYtWlTtZIiIiPQZM3sj33g1y4uIiNQYBXcREZEao+AuIiJSY2rmmruIiPS91tZW1q9fT3Nzc7WTUtOGDh3KmDFjaGpqSjV9RYO7mV0E/AxoAG519x/mfH8u8FPgFOByd5+Z+O6LwPeiwZvcfVol0yplksnA4sUwYQKkLIQiMnCtX7+ekSNHMnbsWMys2smpSe7Oli1bWL9+PePGjUs1T8Wa5c2sAfgl8AngROAKMzsxZ7I3gS8Bd+XMewhwA/BhYApwg5kdXKm0ShmtXQuLFsHChdVOiYj0gebmZkaPHq3AXkFmxujRo0tqHankNfcpwGp3X+vuLcB04NLkBO6+zt2XAJmceS8EnnD3re6+DXgCuKiCaZVyaWsL/1tbq5sOEekzCuyVV2oeVzK4Hw28lRheH40r27xmdpWZLTKzRZs2bepxQkVEZGDavn07v/rVr8qyrLFjx7J58+ayTVdNA7q3vLtPdffJ7j75sMO6PKBHRERqXKnBvS1uXaxxlQzubwPHJIbHROMqPa+IiNSJ6667jjVr1jBx4kSuueYarrnmGk4++WQmTJjAjBkzAHjqqac455xzuOSSSzjxxBNpb2/n29/+NieffDKnnHIKP//5zzuW9/Of/5xJkyYxYcIEVq5cCcCWLVv4+Mc/zkknncSVV16JuwNw/fXX89Of/rRj3u9+97v87Gc/46mnnuK8887jsssu44QTTuALX/hCxzx9pZK95RcC481sHCEwXw58PuW8jwH/O9GJ7uPAd8qfRBERKZerr76axYsXl3WZEydO7BRAc/3whz9k2bJlLF68mPvuu49//dd/5eWXX2bz5s2cfvrpnHvuuQC8+OKLLFu2jHHjxvHrX/+adevWsXjxYhobG9m6dWvH8g499FBefPFFfvWrX3HLLbdw66238oMf/ICzzz6b66+/nocffpjbbrsNgC9/+ct8+tOf5uqrryaTyTB9+nSef/55li5dyksvvcTy5cs56qijOOuss5g/fz5nn312WfOmOxWrubt7G/ANQqB+BbjH3Zeb2Y1mdgmAmZ1uZuuBzwK/MbPl0bxbgf9FOEFYCNwYjRMREcnr6aef5oorrqChoYEjjjiCj3zkIyyM7tyZMmVKx21kc+bM4atf/SqNjaF+e8ghh3Qs49Of/jQAp512GuvWrQNg3rx5/N3f/R0AF198MQcfHOqdY8eOZfTo0bz00ks8/vjjnHrqqYwePbpjfWPGjGHQoEFMnDixY1l9paL3ubv7bGB2zrjrE58XEprc8817O3B7JdMnIiLl010Nu9pGjBiRarohQ4YA0NDQkOr6/JVXXskdd9zBe++9x5e//OUuyyllWeU0oDvUiYhIfRs5ciS7du0C4JxzzmHGjBm0t7ezadMm5s2bx5QpU7rMc8EFF/Cb3/ymI+Amm+XzOffcc7nrrvA4lkceeYRt27Z1fPepT32KRx99lIULF3LhhReWa7N6TcFdREQGrNGjR3PWWWdx8skn8+yzz3LKKafwoQ99iI997GP86Ec/4sgjj+wyz5VXXsmxxx7bMW0cuAu54YYbmDdvHieddBL3338/xx57bMd3gwcP5qMf/Sh/+7d/S0NDQ9m3r6esr3vwVcrkyZNd73PvB1auhHnz4IMfhI98pNqpEZEKe+WVV/iLv/iLaiejajKZDJMmTeLee+9l/PjxFV1Xvrw2sxfcfXLutKq5i4iI9MCKFSs4/vjjOf/88yse2Eult8KJiIj0wIknnsjatWurnYy8VHMXERGpMQruIiIiNUbBXUREpMYouIuIiNQYBXcRERmwqv3K17/8y7/sGH/NNddw0kkncc0117Bp0yY+/OEPc+qpp/LnP/+5LOkrhXrLi4jIgBUH93/8x39MNX1bW1vHM+XL4Zlnnun4PHXqVLZu3UpDQwPTp09nwoQJ3HrrrWVbVylUcxcRkQGrmq98BTjggAMAuOSSS9i9ezennXYaN998M9deey0PPPAAEydOZN++fTz++OOceeaZTJo0ic9+9rPs3r0bgEcffZQTTjiBSZMm8U//9E988pOfLEu+qOYuIiLl8cwzsGVLeZc5ejQkmr5zVfOVr0mzZs3igAMO6Hjl7RFHHMGiRYv4xS9+webNm7npppuYM2cOI0aM4Oabb+bHP/4x1157Lf/wD//A3LlzOf744/nc5z5XtmxTzV1ERGpCX7/yNa3nnnuOFStWcNZZZzFx4kSmTZvGG2+8wcqVKxk3bhzjx4/HzDrWUQ6quYuISHl0U8Outkq98jUNd+eCCy7g7rvv7jQ+ruVXgmruIiIyYFX7la9pnHHGGcyfP5/Vq1cDsGfPHl599VVOOOEE1q1bx5o1awC6BP/eUHAXEZEBq9qvfE3jsMMO44477uCKK67glFNO4cwzz2TlypUMHTqUqVOncvHFFzNp0iQOP/zwkpbbHb3yVcpLr3wVqSv1/srXcnrqqae45ZZbeOihh/J+r1e+ioiI1DF1qBMREekHzjvvPM4777yyLKuiNXczu8jMVpnZajO7Ls/3Q8xsRvT9AjMbG41vMrNpZrbUzF4xs+9UMp0iIiK1pGLB3cwagF8CnwBOBK4wsxNzJvsKsM3djwd+Atwcjf8sMMTdJwCnAV+NA7+IiPQvtdJ3qz8rNY8rWXOfAqx297Xu3gJMBy7NmeZSYFr0eSZwvpkZ4MAIM2sEhgEtwM4KplVERHpg6NChbNmyRQG+gtydLVu2MHTo0NTzVPKa+9HAW4nh9cCHC03j7m1mtgMYTQj0lwLvAsOBb7p7lxsRzewq4Cqg5FsTRESk98aMGcP69evZtGlTtZNS04YOHcqYMWNST99fO9RNAdqBo4CDgT+b2Rx3X5ucyN2nAlMh3ArX56kUEalzTU1NHY91lf6jks3ybwPHJIbHROPyThM1wY8CtgCfBx5191Z33wjMB7rcxyciIiJdVTK4LwTGm9k4MxsMXA7MyplmFvDF6PNlwFwPF27eBD4GYGYjgDOAlRVMq4iISM2oWHB39zbgG8BjwCvAPe6+3MxuNLNLosluA0ab2WrgW0B8u9wvgQPMbDnhJOG37r6kUmkVERGpJRW95u7us4HZOeOuT3xuJtz2ljvf7nzjRUREpDg9flZERKTGKLiLiIjUGAV3ERGRGtNtcDezBjP7Zl8lRkRERHqv2+Du7u3AFX2UFhERESmDNL3l55vZL4AZwJ54pLu/WLFUiYiISI+lCe4To/83JsY50UNmREREpH8pGtzd/aN9kRAREREpj6K95c1slJn92MwWRX//z8xG9UXiREREpHRpboW7HdgF/G30txP4bSUTJSIiIj2X5pr7+939M4nhH5jZ4gqlR0RERHopTc19n5mdHQ+Y2VnAvsolSURERHojTc39a8DvEtfZt5F9TauIiIj0M90GdzNrAP6Lu3/IzA4EcPedfZIyGdjMqp0CEZG61W1wd/f2uEleQV1K4l7tFIiI1K00zfIvmdks4F46P6Hu/oqlSkRERHosTXAfCmyh8xPpHFBwFxER6YfSXHPf4u7f7qP0iIiISC+leSvcWX2UFhERESmDNM3yi3XNXUREZOBI8xCb5DX3v47+Pplm4WZ2kZmtMrPVZnZdnu+HmNmM6PsFZjY28d0pZvasmS03s6VmNjTVFomIiNS5NG+F+/ueLDi6Xv9L4AJgPbDQzGa5+4rEZF8Btrn78WZ2OXAz8DkzawTuJNxj/7KZjQZae5IOERGRepPmrXAfMLMnzWxZNHyKmX0vxbKnAKvdfa27twDTgUtzprkUmBZ9ngmcb2YGfBxY4u4vA7j7luj6v4iIiBSRpln+34DvENWc3X0JcHmK+Y4G3koMr4/G5Z3G3duAHcBo4AOAm9ljZvaimV2bYn0iIiJCug51w939eev8ONG2CqUn1gicDZwO7AWeNLMX3P3J5ERmdhVwFcCxxx5b4SSJiIgMDGlq7pvN7P2EB9dgZpcB76aY723gmMTwmGhc3mmi6+yjCJ331gPz3H2zu+8FZgOTclfg7lPdfbK7Tz7ssMNSJElERKT2pQnuXwd+A5xgZm8DVxPeFFfMQmC8mY0zs8GEpvxZOdPMIvuGucuAue7uwGPABDMbHgX9jwArEBERkaLS9JZfC/yVmY0ABrn7rjQLdvc2M/sGIVA3ALe7+3IzuxFY5O6zgNuA35vZamAr0bV8d99mZj8mnCA4MNvdH+7B9omIiNSdNNfcAXD3PcWn6jLPbEKTenLc9YnPzcBnC8x7J+F2OBERESlBmmZ5kdLpfe4iIlWT5j73IWnGiXSi97mLiFRNmpr7synHiYiISD9Q8Jq7mR1JeMjMMDM7FYjbWQ8EhvdB2mQgU7O8iEjVdNeh7kLgS4T703+cGL8L+J8VTJPUgjTN8vv3w7RpcOaZMGFC5dMkIlInCgZ3d58GTDOzz7j7fX2YJqkXe6IbMFatUnAXESmjNLfCPWRmnwfGJqd39xsrlSipAWqWl3rnDvv2wXBdxZS+l6ZD3QOEt7e1AXsSfyKFqbe81Lvnn4c77wwBXqSPpam5j3H3iyqeEhGRWvLGG+F/czMMG1bdtEjdSVNzf8bMdEFU0ulJjV21fBGRskpTcz8b+JKZvQ7sJ9wS5+5+SkVTJiIiIj2SJrh/ouKpkPqmznciImVVtFne3d8gvHP9Y9HnvWnmkzqlJnYRkapL82z5G4D/AXwnGtWE3tYmIiLSb6WpgX8KuITo9jd3fwcYWclEiYiISM+lCe4t7u6AA5jZiMomSURERHojTXC/x8x+AxxkZv8AzAH+rbLJkgErvuauTnIiIlVTtLe8u99iZhcAO4EPAte7+xMVT5kMbOpYJ9LZ66/D/Pnw+c/DIPVJlspKcyscUTBXQBcR6an582Hv3vDEOj1vXiqsu/e5P+3uZ5vZLqLr7fFXhIfYHFjx1MnA05NmedXyRUTKqmDbkLufHf0f6e4HJv5Gpg3sZnaRma0ys9Vmdl2e74eY2Yzo+wVmNjbn+2PNbLeZfbvE7ZJqU8AWEama7mruh3Q3o7tv7e57M2sAfglcAKwHFprZLHdfkZjsK8A2dz/ezC4HbgY+l/j+x8Aj3W+CDHjqfCciUlbdXXN/gdAcb8CxwLbo80HAm8C4IsueAqx297UAZjad8OrYZHC/FPh+9Hkm8AszM3d3M/sb4HX0etmBSc3yIiJV012z/Dh3fx/h1re/dvdD3X008Eng8RTLPhp4KzG8PhqXdxp3bwN2AKPN7ADCU/F+kHZDpJ9JE7AV1EVEKiLN/RhnuPvseMDdHwH+snJJAkJt/ifuvru7iczsKjNbZGaLNm3aVOEkSSqlBGzdEy/1QCexUgVpboV7x8y+R/Z58l8A3kkx39uEF87ExkTj8k2z3swagVHAFuDDwGVm9iPCZYCMmTW7+y+SM7v7VGAqwOTJk7UH9SdqlhfpTOVc+lCa4H4FcAPw79HwvGhcMQuB8WY2jhDELwc+nzPNLOCLwLPAZcDc6FG358QTmNn3gd25gV36OR3IRALtC1IFaZ5QtxX451IX7O5tZvYN4DGgAbjd3Zeb2Y3AInefBdwG/N7MVgNbCScAMpCpWV6ks9xyrmAvfaBocDezw4BrgZOAofF4d/9YsXmja/Wzc8Zdn/jcDHy2yDK+X2w90g+lCdg6yEk9icu7yr30gTQd6v4ArCTc+vYDYB2hyV2ksJ7U4EVEpCzSBPfR7n4b0Oru/+HuXwaK1tpFilJQl3qQW85V7qUPpOlQ1xr9f9fMLib0lO/26XVSx3py4NI1d6llCu5SBWmC+01mNgr478DPgQOBb1Y0VTLwlXLNXQc7EZGy6ja4R8+HH+/uDxGeHvfRPkmVDHwK2CKBau5SBd1ec3f3dtLd0y4S6FY4ke4puEsfSNMsP9/MfgHMIPESF3d/sWKpkoFPT6gTCVS+pQrSBPeJ0f8bE+Mc9ZiX7ujFMVLvCpVvlXvpA2meUKfr7FJZapaXWqZr7lIFaZ5Q9608o3cAL7j74rKnSAa2Uq6jq7e81AOVb6mCNA+xmQx8jfDu9aOBrwIXAf9mZtdWMG0ykOmAJpKf9g3pA2muuY8BJsXvVjezG4CHgXOBF4AfVS55UtPUW17qkYK79IE0NffDgf2J4VbgCHfflzNepGcBWwc7qWUq31IFaWrufwAWmNkD0fBfA3eZ2QhgRcVSJgObesuLBOpQJ1WQprf8/zKzR4CzolFfc/dF0ecvVCxlIiIi0iNpau5EwXxR0QlFYqX0ltc1d6llqrlLFaS55i5SOr3PXSSdlSvhzjsrv54dO2Dr1sqvp79yhz/9CTZtqnZK+oSCu5SXgnr/89ZbsHFjtVNRv4rV3OfNg717K78/zJgBM2dWdh392d698Npr8Nhj1U5Jn1Bwl8pQs3xWa2t11//II/DHP1Y3DcXs3BlOQuqBHksrfUDBXcqrJ0+dq+WD2iuvwG9/G4KXFDZ9ejgJ6e+am0ufJ235ruX9QPqcgrtIJb3+evi/Y0d10yG99+qr8LvfwebNvVuOau7SByoa3M3sIjNbZWarzey6PN8PMbMZ0fcLzGxsNP4CM3vBzJZG//UGulpUL83yUhvefjv837attPnS9pbPZEpPk0gBFQvuZtYA/BL4BHAicIWZnZgz2VeAbe5+PPAT4OZo/Gbgr919AvBF4PeVSqdUSCkPsamHGotOYKSYetgPqqnOKhOVrLlPAVa7+1p3bwGmA5fmTHMpMC36PBM438zM3V9y93ei8cuBYWY2pIJpHfg2b4Zbbw09QgWmTu0fvWJ1wJa0NXeVlcqqs/ytZHA/Gkh2f10fjcs7jbu3EV4lOzpnms8AL7p7l+fYm9lVZrbIzBZtqpN7Fwtatiw061W7x3F/2oHeeKPaKZB6VqhlSsG9OuqppZB+3qHOzE4iNNV/Nd/37j7V3Se7++TDDjusbxMnvVcPO1mdNQVKL/Tl/tDWVhvX+PfsCduSRjK49+SuhwGmksH9beCYxPCYaFzeacysERgFbImGxwD/DvxXd19TwXRKJejFMZ0puKczEAJOqeW2P9bcb78dnnii79ZXKX/4A8yenW7aOH/37Qt3PaQ9KRigKhncFwLjzWycmQ0GLgdm5Uwzi9BhDuAyYK67u5kdRHhn/HXuPr+CaZRyG6gBe+/ecJ1+w4byLneg5odUTn8I7lA7l63eey/ddLknju3t5U9LP1Kx4B5dQ/8G8BjwCnCPuy83sxvN7JJostuA0Wa2GvgWEN8u9w3geOB6M1sc/R1eqbRKlfSnwBff5rR8eXXTUe8GQs29VHqITf9QZ/mb6q1wPeXus4HZOeOuT3xuBj6bZ76bgJsqmTapsDrbkQqqs048vVYP+dRfau615PXX4dBDYeTIwtPU2dv5+nWHun6rrS1ct5HeqfGdq5N62tbeGAj5VGr/CT3EpjKS+fXEE3Dffd1Pn5vvNZ7fCu498cc/wu/1XJ28enJwrofOZgMhaPUH9ZxPfbXttZLHudvR0lLa9LWSDwUouPdEPb8TuZzqocm6HraxnGoxn/pbb/laqbEWyq977oEHHyw+fS2WtQQFdwn3ipZrh9+1qzzL6Y3e7LSVakWo1gG1Ugewp5+GF18s/3JrJfD0hGruwe7d6SpQhcrK9u3w7rtdx6tZXupKe3u4V3TevPIsb+3a8L+a97n3p4NXnJZMpjotPuXOi0wGNm6EFStg0aLyLhv6129XLqq5l+auu2DmzOLTlbodudP393zoJQX3ehff6xm/mrQ3+suBuSfpqHTan38+HLC2b6/senKVe7ueeSb0OamU/lKGutPbh9iUa7k9NRDyOI1KPUyoRii417tynr2Wuqx6qLnHdu4M//v6xT7lzovevsu8mP742+XqbRrdw+Wr3Jacvtr2Wnl4S7H8yj0eqVleUhsIB6Ji4gJejmvNyZ2lmnnTn3baapeR/pQXaVQ7v9IoR3C/++6uTc/J3+rBByvXQtIXZeL118MTHyvZB6fYduzPedeYau6SWn88cJZaYMu5DaXWCOqp5l4tafLiueeyfSWqrT/uU9A5UFSqWT7p3XdD34ZK6Iua+2uvhf/lbunJZGD+/NDprlhZyb01rho19z17it+iVyEK7r3RHw9E+dJ0772hA1Ta6cu57mroSTr66kTDPTTHTp3aN6/nTbNdS5bAnDmVT0safXFitnNnyP+0zyRfswamTYP4tdJp07hnT/7pq/0Qm75sli/X3Sdr14b9ZvPm8IjouXOL/w6trZ2Hq1Fz/8Mfij9cp0IU3HujvwSzpHwFdtu2cOtSPvE2lKOgD+Rr7n35atY4qKxbV3iaNWvKU3NL5sXtt6ebZ9++EPy6S1+l9MUBN36PwKuvljZ93BkyTTl/773CQXSg9pZft678L1ZKa86ccBmjoSEMNzcXz69K1Nx37Ch9nirdHqzg3hsDIbgXS2Olau4D7Va4vvot067nySfLc801mRfJV1zu31+4c9+WLeF/NV6i09MysWFD5WqkPanxbdtW2vSlTJdr0aLO6yump/n0+OPwwAM9m7dc4jxqaSm9Wb5YB7tiNmyAGTMKt4L2MwruvdEfr+2W+lpDNcv3fJ40ujvZ6otWgkLbdeed4a+/6W6f2rMn/4F1584QdJ55Jt06cvN9x478Dz3pSRpjg7o5tOZ29CplublaWsLDhB56KP081ewtH981kka+/IjTnia4F2uWL3Wfj2vtleoLUWYK7r3RX4JZUm6BrsfgXugguXcvPPVU5xpssXnKnZb29r49KSy0rrhc9GWLxdy5xe/z7y49jzwSLi/lvrSpuTn8L7XzVpw3M2bkf1xp7nSFhvNJnkDkTv/ss9nPvb3DJJ6/2H6eXHacX0m7dmUvP1TKhg0wfXr6mm+8bck8irezrS1/fiXHlTu490VHyjJScO+NTAbefDPUgCp5Nrx7d/4dMp8XXug8XOngvnt3Nlj2l2b5Qtv0/PPhOuvKlSGfkusfSJ2Zli3L1v7cu19msTzuq1rc1q2wenUI8N3JF0hXrw6/T7wPVPs+7d7W3JOS29KT/SHfiWo+yWXn6zx5773w8MOlr7+QuL9Gcr+KT+zS1nzjeZPbWCy/kusrds09Hn777Z5dSy+myuVUwb03MplwFr53bwhyaeQejPftK14I7rqrtLfQFdoZ8ultAbzrLpg9O3weKB3qnnkmBPc33siOK5b2vXvTPz7WPZxI7N3bNX8zmeKd93LT8swz4ZWWEGo/zzwDf/5zGJ47F267LTvt/v1h/vhAWqy2ku/3783vEl+vzxVva25tKldu+tasCdu4ZEl2XG5Aq2QNqaWla6/6NOU8bdAuV3DvroVm9+7iy057kpBGMr960zKRr1Ui+Tnf79CT4P7ww6H1ZubM8r4vocotmQruvdGTgvvnP3c+GP/+96G5sZhSdozkATRtzb21tfiBt9C88c5caF19/cjV3NaLWG4wLeX3e/DBsPPnq03kevddWLw4NCHn5kmak6ncaZYtCw8FyWSy96PHTdNr1mSnW7o03LI1d254M1a+g3qa4B6PK6XMbd8e0nnfffnvmY/zq1gZe/DBzh394haK5LXa3LwvNTDF25W7fVu3dl3WI490vU7sHloRuvstexLc0zaxF5o/n//4j3ACnnaZhfKylECVXMamTSFfV60q/f0VxWruxYJ7qc3yW7eme19C2r4yuXm+Ywe88066ectAwb2YNWsKP+AjWThyd4p9+/LfSrRyZfjf2pqdP/7Bm5vDAbonktfLkmespTTL33tvaevM3eZ8O9tbb4VAs3p1acvuqd27C/9e3e2UyVsC29pg/frO38fNdhs2hPy9/XZ46aXO07S2hgNE8qCUL4962pqyeHHh8tHamn0/QPIEIPc3yU1Payv/8i//wuNxy0C+adyLH9zvuSfboS2+9t3Skj0JiZeZJhDHwXTp0vDAknhZ8e+Xu4xSW58KTT9zZgiGSflu/XKH3/0OHnus8DrSnnB0d3IFxfM933peey2b73FZSHtZr9DJV3fb0913S5eGfE3edpg2OBaruRdrlk8T3Ivl7zvvZDtblloTz/1tn346dHxM+3yFXlJwL+bJJws/4OP++7OBNFmQ9u8PNfLHH8+Ob27ueitSbq/ZefNCM3/8sIy0du7MPjADOgf3fDteslNXssDu3h2W9fLLYTi+Nl1I7s6TryYcN2Xn6+xUiebU7g40uQeVOXPCwXvTpuxB3D3shLNnd25xGDw4/N+3L3vgXLWq8/IefbRz7T5fetrbszt9oYNcoW1I1ir27+988Mh3UE5eAkiuP2n/fpYtX87M5KNQc6f505/g1lvzpymfePvvvjt7Oam74J7bwSpOc3winDtf7rbG3/W0RpWUplNZvL7cE8C06yg0XdwDvJTgnruePXvC7zVrVjjZi6/99yS4ZzLheQcvv1x4e958M5zoJvfvNEEwebx68cX8rXvJznP5lp3mmntzc/bEPN+tcMVOwh56KLQm7dpV+nX53DyLjxtHHlnacnpIwb0U77zTtaduPJzcKZIHqziA/+53nTusNDd3Du5tbfmXldsBqa0tHORzTxSSkkEnt4C5h8sCzz0XhnML/KOPwoIFoWl03ryeB/c0cnfOlpaQd6XclpQrNy+SJzr5Dv7PPht24Pigvn599p7h5LLieffvL/w4yTjdySeTdRfcC+VX8jcr9NCQrVs7X97JF9xbWvLncUJr8iCbb/2QbXUp1gyaOy6Zf3H68r12M/cBS4W2pVDNvdRm+d72M8k9BhRbR3cnscn70+fNCyecpfSZyd32+PfdsSP004i/T5Nm6Jz3ccBdsqRwHsf9VpLlNN/+kcyDjRvDk9tWrw7HwUWL8t/Kl69ZPvk5btVJSp5wtbSEFqUZM7qmIV5+oZaK7duzfYkgnKgm+320txfvGJj727W0wPjx3c9TRhUN7mZ2kZmtMrPVZnZdnu+HmNmM6PsFZjY28d13ovGrzOzCSqYzld27QwH805/yf58sJPFTlCAU3riQJneA3Jp78sw6eWtPbnP20qXhTHfFimzzZW5AS55c5Baw+Hpm3Lyb+3287uTBILezYLyT5DZLdhfc3cMyu2uef+yxcKDPvS1pyZLsjrZxY+FbnpYs6ZqmZODId5DduLHzb9fcnD3gJn+TuAa0ZUvxA2X8uxRqco0PUIUO3MnxaR8aku+Ami+4//u/dzoo7ch3gCp0XTDfYzTnzOl8axcUryGtXZsdztdJMd6W5MlY8jcqFNw3bux6SWbOnND6ltTb2wCT/Rxefz10nsy1bFn2c3d9TnJbBdet69zn4OWXu382eby9ra3huFCoht5dmU02EyfzOd4Phg0rfpKR/K2KnWzFy507N1R6IGzzu+/mv4Ml2Ukz2XqVryadvKyyaVM2P9zzd6gr1OI1f373LTPPPhseMvXee4UfI52nlYyhQwsvs8waK7VgM2sAfglcAKwHFprZLHdPtsF9Bdjm7seb2eXAzcDnzOxE4HLgJOAoYI6ZfcDd++7egr17O19Tveuu8L9Q00xra/ihd+/uGizyPQmsuRmamjoPx4Uvnn7Bgs7z7NiR3dFffTXUvi+8sPvr9LkFLPdRiIUOcMn3u2/dCgccENL34INhO6+8smteJA9C8bbEebF0aTadRxwBI0dmp9m2LRzEkjX2114LZ7nvvZdtZdiwIRvsjjoKzjsvpCsWT5cU73gbN6Z/gEZ8cHruuXCiNmZMNh9XrAhNkbmS+Rhf1sh3kHv55ewOvmoVHHooDB8e8mDSpMLzFZPvIPXaa3DMMZ3HtbV1evJdp+C+c2fYzmSATJ6M5XsKWr5+JW++2TmP9u/vvE1z5sAHPwhnnpn/N8kXzJIH+ELN8vGyr7oq7E9LlmS35fzzs9PEv2V3nfv27StcM0v+1nFfhUmTQjB58EE47bTO27BqFZx+euGH1+RK7leLF4e/K6/MnmBmMtnPyWvZixbBlCn5l1nohG3DhtCEH4vzpK0tuz8OGdK19tzQUDig5/v90lxnfvNNOOyw7HCcz8kaeikvYUn+voU6l+YrAwsWdF82krX2OO8+85kQ8CdODMe2+K6V3HUNGZI+/b1kXqHbSMzsTOD77n5hNPwdAHf/P4lpHoumedbMGoH3gMOA65LTJqcrtL7Jkyf7ojQ9HVPYtX07S6++Ot9GdR52J86/9sZGGqICbvF0ZrQOG8begw5iVE5Tc8uIETQfcAAHRrX5TEMDg6Jr4ZmGBtoGD2ZwnrPttiFDaNy/v+j1xW1jxjBy82asvZ2GqKDuPeQQhu7a1TG888gjOTBl546NH/gAIzdsYFh04Nly3HEckriVbNsxx3BwTtBrHTaMptxtiNK9/eijOejtt7ttstx55JEd+QN0mjbO9+YDD2Tn4YczuLmZg6IzbSvjk98yjY1sev/7OfzVV7GctGYaGth2zDE07d3LyC1bsJwTJR80qMu47uwfMYIhe/awZ/RoRhS6payA9qamjt81m4Di+/amTZt4OGoV+fsvfanL9+U6OuRNX0q5v2bzyJG0DRvGARs3su+ggxjU3s6QxEnr5nHjGLxnDyMTwXnL2LFkGhsZtmMHTfv2MWT3bjBj+1FHhXKYsPvwwzvNm08lH0/SMnw4g/NUCLYedxwAB61fz56DD6Zl+HBGJ2/njPWy/G859tiuy82zzJbhwxkcXdbZN2oUzQceyLAdOxiaPGHLmS851CUPc8rrluOOo72picNXr+5YTtot6zieFpEvr3tTVpOaR41iaHS83Dp2LIesW8fG8eM552tf45BDDun18mNm9oK7T84dX7GaO3A0kGyvWA98uNA07t5mZjuA0dH453LmPTp3BWZ2FXAVwLHHHlu2hG9Zu5ap06YxOGf8GuB9pCtgG4GDgaYi020FRgENwC5gGXBmMi2EDMlnFzAyMbwSaAFOiYbbo+Xm2gsMI/2O0h9sBIZHf1sIZ4DFrAHen2c5h/cwDXuAESXO00Zld7Le2EYoe/G1ud/ecQcZKnutLl+ZrPQ682khBJdi9aj3gJ50f8oA04HPJ8bl7q87CPkP0AwMTXxeCnwQOLAH6861GTiUwseDStoNvAmcmDN+I2HbhpDuOPQGcFzKdW4HFgPnUdr+3gKdjvm7CfvvQcCrhP34fXnmS7uOvcAfgfkXX1zW4F5Ifz3upOLuU4GpEGru5Vru0RMm8P01a6CtDctkcDO8sTE0he3fD42NeHs71tREQ0sLNDZimQyZpibcHWtpIdPU1Pl+4cGDQ9NV3Jy2bx8G+NChoYkr0dRmra3hszs0NoIZtm8f3tAQ1mOGuZNpamJQSws+ZEiYNp4/auL3YcOyy21txdvbwxlwY2P28Y1ROi2TwZuaQtPR4EQRj24Xsag5zAcPhkGDsLY2zCwMR9vsgwZhcXpbWsK2tbeH9be342bhPzAokwnpJtq5zTo/3GXvXnzQIHzwYKy9PaQZQk24oQHa2xnU2oo3NWHRf8zCOiA03yWbExsasObmkKbW1s79IuJ8jq4NW3t7drva2sI2mIX8jK/TRdvIoEGhfAweHNYZX2rJZLLpGjQo5GNDA7Z/f0gDZJs34/ni36qtLVvjcQ/bP2hQxzaSyeBELRSZTBgXL6uhIaQvmrajfMTLilp93Cz8zu3tjBgyhL379tE+aFDH8oGw3Q0NWJxX+/eH9MXbGJfT+O6LwYOzvb2jMuJDh0JLC+aOx+mLy2q8f8R5EJVHa2kJ25rJhP0uzvPGxrDd8V0ADQ3ZfaqxsVM5tUwmW/4zmdDqEpVxS+4z7e0d+0DHfh3/RtE6PPoNrLU1zBf9ThbvT1G5NfeQ7ihvvhtPG6fRLLvdcbO2e7asxtsZ/87x/uce0pXYxwEGtbfj0bZ35FlbWxgXHYsGtbSEfai5OeR/U1N2PfFvEZXlzJAhYV9sa8OiskFjI+zfjzc24mZhf4j28bicWrycuCxG+7BnMh37bbzfd5SrKJ86WtmivIz3b2tp6Vh+ZvDg8NtkMtm8istclH7cs8ewOI/itLa0ZPM7kT6iY3cyPfF+S1x+EvtBx+/R3h6OEY2N4f+QIdljQltbKLPQ8VtblNfe1MSNgwZx9NFd6qkVUcng/jaQvOA3JhqXb5r1UbP8KELFLM28FdPU1MTY9+U7RxMREen/KtkathAYb2bjzGwwoYPcrJxpZgFfjD5fBsz1cDF1FnB51Jt+HDAeyNMlVURERHJVrOYeXUP/BvAY4VLP7e6+3MxuBBa5+yzgNuD3ZraacPn58mje5WZ2D7CCcNnj633aU15ERGQAq1hv+b5Wzt7yIiIiA0Gh3vJ6Qp2IiEiNUXAXERGpMTXTLG9mmwi3Q5bToYTbRKV0yrueU971jvKv55R3PVetvDvO3bs8+qNmgnslmNmifNcypDjlXc8p73pH+ddzyrue6295p2Z5ERGRGqPgLiIiUmMU3Ls3tdoJGMCUdz2nvOsd5V/PKe96rl/lna65i4iI1BjV3EVERGqMgnseZnaRma0ys9Vmdl2109PfmNkxZvYnM1thZsvN7J+j8YeY2RNm9lr0/+BovJnZv0T5ucTMJlV3C6rPzBrM7CUzeygaHmdmC6I8mhG9j4Ho/QozovELzGxsVRPeD5jZQWY208xWmtkrZnamyl46ZvbNaJ9dZmZ3m9lQlb3CzOx2M9toZssS40oua2b2xWj618zsi/nWVW4K7jnMrAH4JfAJwmuIrzCz3NcR17s24L+7+4nAGcDXozy6DnjS3ccDT0bDEPJyfPR3FfDrvk9yv/PPwCuJ4ZuBn7j78YRXrX8lGv8VYFs0/ifRdPXuZ8Cj7n4C8CFCPqrsFWFmRwP/BEx295MJ7/y4HJW97twBXJQzrqSyZmaHADcAHwamADfEJwSVpODe1RRgtbuvdfcWYDpwaZXT1K+4+7vu/mL0eRfh4Ho0IZ+mRZNNA/4m+nwp8DsPngMOMrP/1Lep7j/MbAxwMXBrNGzAx4CZ0SS5eRfn6UzgfOt4CXb9MbNRwLmEl07h7i3uvh2VvbQagWHRK7aHA++isleQu88jvNQsqdSydiHwhLtvdfdtwBN0PWEoOwX3ro4G3koMr4/GSR5RU92pwALgCHd/N/rqPeCI6LPytLOfAtcCmWh4NLDd3dui4WT+dORd9P2OaPp6NQ7YBPw2uqxxq5mNQGWvKHd/G7gFeJMQ1HcAL6CyV6pSy1pVyqCCu/SYmR0A3Adc7e47k995uA1Dt2LkMLNPAhvd/YVqp2WAagQmAb9291OBPWSbRQGVvUKipuBLCSdIRwEj6IMaZC3rz2VNwb2rt4FjEsNjonGSYGZNhMD+B3e/Pxq9IW7yjP5vjMYrT7POAi4xs3WESz4fI1xDPihqKoXO+dORd9H3o4AtfZngfmY9sN7dF0TDMwnBXmWvuL8CXnf3Te7eCtxPKI8qe6UptaxVpQwquHe1EBgf9SAdTOhwMqvKaepXoututwGvuPuPE1/NAuKeoF8EHkiM/69Rb9IzgB2JZq264u7fcfcx7j6WULbmuvsXgD8Bl0WT5eZdnKeXRdP3y5pCX3D394C3zOyD0ajzgRWo7KXxJnCGmQ2P9uE471T2SlNqWXsM+LiZHRy1nnw8GldZ7q6/nD/gPwOvAmuA71Y7Pf3tDzib0BS1BFgc/f1nwvW4J4HXgDnAIdH0RrgDYQ2wlNBbt+rbUe0/4Dzgoejz+4DngdXAvcCQaPzQaHh19P37qp3uav8BE4FFUfn7I3Cwyl7qvPsBsBJYBvweGKKy121+3U3on9BKaDX6Sk/KGvDlKB9XA3/fF2nXE+pERERqjJrlRUREaoyCu4iISI1RcBcREakxCu4iIiI1RsFdRESkxii4i0jFmNlTZja52ukQqTcK7iIiIjVGwV2kzpjZCDN72Mxejt7r/Tkzu97MFkbDU+O3f0U175+Y2SIL704/3czuj95LfVM0zVgL71b/QzTNTDMbnme9HzezZ83sRTO7N3o3AWb2QzNbEb0D+5a+zQ2R2qTgLlJ/LgLecfcPeXiv96PAL9z99Gh4GPDJxPQt7j4Z+FfCoza/DpwMfMnM4reEfRD4lbv/BbAT+MfkCs3sUOB7wF+5+yTCE+a+Fc3/KeAkdz8FuKkymyxSXxTcRerPUuACM7vZzM5x9x3AR81sgZktJbzM5qTE9LMS8y1393fdfT+wluwLMd5y9/nR5zsJjyhOOgM4EZhvZosJz+Q+jvAa0WbgNjP7NLC3nBsqUq8ai08iIrXE3V81s0mE9wHcZGZPEmrjk939LTP7PuG54rH90f9M4nM8HB9Dcp9jnTtswBPufkVuesxsCuElJpcB3yCcXIhIL6jmLlJnzOwoYK+73wn8X8IrUwE2R9fBLys4c2HHmtmZ0efPA0/nfP8ccJaZHR+lYYSZfSBa3yh3nw18E/hQD9YtIjlUcxepPxOA/2tmGcLbrv4b8DeEN4W9R3jtcalWAV83s9sJrxH9dfJLd99kZl8C7jazIdHo7wG7gAfMbCihdv+tHqxbRHLorXAi0itmNpbw6tqTq50WEQnULC8iIlJjVHMXERGpMaq5i4iI1BgFdxERkRqj4C4iIlJjFNxFRERqjIK7iIhIjVFwFxERqTH/HyuKwvgSV4dWAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize=(8,3))\n",
    "ax = plt.subplot(111)\n",
    "ax.plot((grad1-grad_true).abs().sum(1), c='black')\n",
    "ax.plot((grad2-grad_true).abs().sum(1), c='red', alpha=0.4)\n",
    "ax.set_ylabel('gradient error')\n",
    "ax.set_xlabel('samples')\n",
    "plt.legend(['torchdyn', 'torchdiffeq'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Interpolated Adjoint"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.03173661231994629\n",
      "0.019273042678833008\n"
     ]
    }
   ],
   "source": [
    "t0 = time.time()\n",
    "prob = ODEProblem(f, sensitivity='interpolated_adjoint', interpolator='4th', solver='dopri5', atol=1e-4, rtol=1e-4)\n",
    "t_eval, sol_torchdyn = prob.odeint(x, t_span)\n",
    "t_end1 = time.time() - t0\n",
    "print(t_end1)\n",
    "\n",
    "t0 = time.time()\n",
    "sol_torchdiffeq = torchdiffeq.odeint_adjoint(f, x, t_span, method='dopri5', atol=1e-4, rtol=1e-4)\n",
    "t_end2 = time.time() - t0\n",
    "print(t_end2)\n",
    "\n",
    "true_sol = torchdiffeq.odeint_adjoint(f, x, t_span, method='dopri5', atol=1e-9, rtol=1e-9)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.8546793460845947\n",
      "0.844252347946167\n"
     ]
    }
   ],
   "source": [
    "t0 = time.time()\n",
    "grad1 = torch.autograd.grad(sol_torchdyn[-1].sum(), x)[0]\n",
    "t_end1 = time.time() - t0\n",
    "print(t_end1)\n",
    "\n",
    "t0 = time.time()\n",
    "grad2 = torch.autograd.grad(sol_torchdiffeq[-1].sum(), x)[0]\n",
    "t_end2 = time.time() - t0\n",
    "print(t_end2)\n",
    "\n",
    "grad_true = torch.autograd.grad(true_sol[-1].sum(), x)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f7829a303d0>"
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAf4AAADQCAYAAADmvgPXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAA8y0lEQVR4nO3dd3gc1dnw4d+jYsmWuzDFGHcHYxvjIlowDuWlt0BMcWihBEhMCAkxgZCYQMgb4KWEEEz5gOCQBEwCAQcMNsY4xgYbV9yL3GTJTcVWb7t7vj9mVhpLq9WstFX73Nela7WzM7Nnz5w5z5wzZ2bEGINSSimlkkNKrBOglFJKqejRwK+UUkolEQ38SimlVBLRwK+UUkolEQ38SimlVBLRwK+UUkolkbRYJyAajjjiCDNw4MBYJ0MppZSKihUrVhQZY/oE+iwpAv/AgQNZvnx5rJOhlFJKRYWI7GrpM+3qV0oppZKIBn6llFIqiWjgV0oppZJIUpzjD6S+vp78/HxqampinZQOKzMzk379+pGenh7rpCillLIlbeDPz8+nW7duDBw4EBGJdXI6HGMMxcXF5OfnM2jQoFgnRymViLZsgW3b4KKLYp2SDiVpA39NTY0G/QgSEbKzsyksLIx1UpRSier737deNfCHVVKf49egH1mav0opFX+SOvDH0qFDh5g+fXpY1jVw4ECKiorCNp9SSqmOSwN/jIQa+D0eTwRTo5RSKllo4I+RBx54gG3btjFmzBimTp3K1KlTGTVqFCeeeCIzZ84EYMGCBZx55plcfvnljBgxAq/Xyy9+8QtGjRrF6NGjef755xvW9/zzzzNu3DhOPPFENm3aBEBxcTHnn38+I0eO5Pbbb8cYA8C0adP44x//2LDsQw89xHPPPceCBQs466yzmDRpEsOHD+f6669vWEYppVTHkLSD+5zuvfdeVq9eHdZ1jhkz5rDg2tTjjz/OunXrWL16Ne+++y4vvfQS33zzDUVFRZx88slMnDgRgJUrV7Ju3ToGDRrEiy++yM6dO1m9ejVpaWmUlJQ0rO+II45g5cqVTJ8+naeeeopXX32VRx55hAkTJjBt2jQ++ugjXnvtNQBuvfVWrrrqKu699158Ph9vv/02X3/9NWvXrmXVqlWsX7+evn37csYZZ7B48WImTJgQ1rxRSikVO9rijwOLFi1i8uTJpKamctRRR/Gd73yHZcuWAXDKKac0XA43b9487rzzTtLSrOO13r17N6zjqquuAmD8+PHs3LkTgIULF3LDDTcAcMkll9CrVy/AOtefnZ3NqlWrmDt3LmPHjiU7O7vh+/r160dKSgpjxoxpWJdSSqmOQVv8ELRlHmtZWVmu5svIyAAgNTXV1XiA22+/nTfeeIN9+/Zx6623NltPKOtSSimVOLTFHyPdunWjvLwcgDPPPJOZM2fi9XopLCxk4cKFnHLKKc2WOe+883j55ZcbgrGzqz+QiRMn8o9//AOAjz/+mIMHDzZ8duWVV/LJJ5+wbNkyLrjggnD9LKWUUnFOA3+MZGdnc8YZZzBq1Ci++uorRo8ezUknncQ555zDk08+ydFHH91smdtvv53+/fs3zOsP6i15+OGHWbhwISNHjuS9996jf//+DZ916tSJs88+m2uuuYbU1NSw/z6llFLxSZJh1HZOTo5Zvnz5YdM2btzICSecEKMUxZ7P52PcuHH885//ZNiwYRH7nmTPZ6VUO+TkWK9N6m/VOhFZYYzJCfSZtviT0IYNGxg6dCjnnntuRIO+Ukqp+KOD+5LQiBEj2L59e6yToZRSKga0xa+UUkolEQ38SqmOq7QUFi6MdSqUiiva1a+U6rh+9jNYswbmz4fu3WOdGqXigrb4lVIdV36+9ao3olKqgQb+GIn1Y3m//e1vN0yfOnUqI0eOZOrUqRQWFnLqqacyduxYvvjii7CkTymlVPzQrv4Y8Qf+H//4x67m93g8DffoD4cvv/yy4f9XXnmFkpISUlNTefvttznxxBN59dVXw/ZdSiml4oe2+GMklo/lBejatSsAl19+ORUVFYwfP54nnniC+++/nw8++IAxY8ZQXV3N3LlzOf300xk3bhxXX301FRUVAHzyyScMHz6ccePGcc8993DppZdGK+uUUkq1g7b4AZ5+GjZvDu86jz8e7ruvxY9j+Vhep1mzZtG1a9eGxxIfddRRLF++nD//+c8UFRXx2GOPMW/ePLKysnjiiSd45plnuP/++/nhD3/I/PnzGTp0KNdee214804ppVTEaIs/DkT7sbxuLVmyhA0bNnDGGWcwZswYZsyYwa5du9i0aRODBg1i2LBhiEjDdyillIp/2uKHoC3zWIvUY3ndMMZw3nnn8dZbbx023d87oJRSKvFEtMUvIheKyGYRyRWRBwJ8niEiM+3Pl4rIQMdnD9rTN4vIBfa0TBH5WkS+EZH1IvJIJNMfSbF+LK8bp512GosXLyY3NxeAyspKtmzZwvDhw9m5cyfbtm0DaHZgoJRSKn5FLPCLSCrwAnARMAKYLCIjmsx2G3DQGDMUeBZ4wl52BHAdMBK4EJhur68WOMcYcxIwBrhQRE6L1G+IpFg/lteNPn368MYbbzB58mRGjx7N6aefzqZNm8jMzOSVV17hkksuYdy4cRx55JEhrVcppVTsROyxvCJyOvBbY4y/tf4ggDHmD4555tjzfCUiacA+oA/wgHNe53yOZbsAi4AfGWOWBkuLPpY3shYsWMBTTz3Fhx9+2OwzzWcVU+efDyUlMHcuOMbEqAShj+Vts1g9lvdYYLfjfb49LeA8xhgPUApkB1tWRFJFZDVwAPi0taCvlEpiEWrYKJXIEm5UvzHGa4wZA/QDThGRUYHmE5E7RGS5iCwvLCyMahqTzVlnnRWwta+UUir+RDLwFwDHOd73s6cFnMfu6u8BFLtZ1hhzCPgcawxAM8aYV4wxOcaYnD59+rT9VyilEpdIrFOgVNyJZOBfBgwTkUEi0glrsN6sJvPMAm62/58EzDfWoINZwHX2qP9BwDDgaxHpIyI9AUSkM3AesKmtCYzU+AZl0fxVSqn4E7Hr+I0xHhG5G5gDpAKvG2PWi8ijwHJjzCzgNeBNEckFSrAODrDnewfYAHiAKcYYr4gcA8ywR/inAO8YY9rUx5yZmUlxcTHZ2dmItgrCzhhDcXExmZmZsU6KSmZ68KlUMxG9gY8xZjYwu8m0aY7/a4CrW1j298Dvm0xbA4wNR9r69etHfn4+ev4/cjIzM+nXr1+sk6GUHgAo5ZC0d+5LT09vuBWuUqqD8vfmaeBXqkHCjepXSqmQaeBPbLr9wkoDv1Kq49PAoVQDDfxKKaXimx64hZUGfqWUUiqJaOBXSnV82mJUqoEGfqVUx6eBP7Hp9gsrDfxKqY5LA0bHoNsxrIIGfvtJeD+LVmKUUioiNHAo1SBo4DfGeIHJUUqLUkqFl97Ap2PQ7RdWbu7ct1hE/gzMBCr9E40xKyOWKqWUUspPA39YuQn8Y+zXRx3TDHBO2FOjlFLhpAFDqWZaDfzGmLOjkRCllIoYPQBQqkGro/pFpIeIPCMiy+2/p0WkRzQSp5RS7aLn+DsG3X5h5eZyvteBcuAa+68M+EskE6WUUmGlgSOx6fYLKzfn+IcYY77neP+IiKyOUHqUUip8NGAo1YybFn+1iEzwvxGRM4DqyCVJKaXCTA8AEptuv7By0+K/C/ir47z+QeDmyCVJKaXCxH+OXyU2DfxhFTTwi0gqcKMx5iQR6Q5gjCmLSsqUUipcNHAo1SBo4DfGeP3d/BrwlVIJRwO+Us246epfJSKzgH9y+J373otYqpRSKpz0ACCx6fYLKzeBPxMo5vA79RlAA79SKjFo4FCqgZtz/MXGmF9EKT1KKRU+egOfjkG3X1i5eTrfGVFKi1JKKdWcBv6wctPVv1rP8SulEpIGDKWa0XP8SqmOTw8AEptuv7By83S+W6KREKWUCjs9x69UM26ezvctEflMRNbZ70eLyK8jn7QOprAQNm2KdSqUSk4a+JVq4OZe/f8PeBCoBzDGrAGui2SiOqQrroAbboh1KpRKLhrwOwbdjmHlJvB3McZ83WSaJxKJ6dDq6mKdAqWSVywDR1mZ7v/tpYE/rNwE/iIRGYI1oA8RmQTsjWiqlFIqHOLhIT3nnANTpsQ6FYlNA39YuRnVPwV4BRguIgXADuD6iKZKKaXCKdaBY9Wq2H6/Ug6ttviNMduNMf8D9AGGG2MmGGN2uVm5iFwoIptFJFdEHgjweYaIzLQ/XyoiAx2fPWhP3ywiF9jTjhORz0Vkg4isF5Gfuv6lSqnkE+uAr8JDt2NYuenqB8AYU2mMKXc7v3273xeAi4ARwGQRGdFkttuAg8aYocCzwBP2siOwBhCOBC4Eptvr8wD3GWNGAKcBUwKsUymlDqeBQ6kGrgN/G5wC5No9BnXA28AVTea5Aphh//8v4FwREXv628aYWmPMDiAXOMUYs9cYsxLAPgjZCBwbwd+glEpkeh2/Us24uY4/w820AI4Fdjve59M8SDfMY4zxAKVAtptl7dMCY4GlLaT7DhFZLiLLCwsLXSRXKdVhaeBPbLr9wspNi/8rl9OiRkS6Au8C9xpjygLNY4x5xRiTY4zJ6dOnT3QTqJSKDxowOgbdjmHV4qh+ETkaq5XdWUTGAv7rYroDXVysuwA4zvG+nz0t0Dz5IpIG9MB6LkCLy4pIOlbQ/7s+KEgppZQKTbDL+S4AfoAVdJ9xTC8HfuVi3cuAYSIyCCtoXwd8v8k8s4CbsXoQJgHzjTHGfhrgP0TkGaAvMAz42j7//xqw0RjzDEop5Ya2GBObbr+wajHwG2NmADNE5HvGmHdDXbExxiMidwNzgFTgdWPMehF5FFhujJmFFcTfFJFcoAT7VsD2fO8AG7BG8k8xxnhFZAJwI7BWRFbbX/UrY8zsUNOnlEoCOrivY9DtF1ZubuDzoYh8HxjonN8Y82hrC9oBeXaTadMc/9cAV7ew7O+B3zeZtojGUw5KKRWcBgylmnET+D/AGm2/AqiNbHKUUioC9ABAqQZuAn8/Y8yFEU+JUkpFigb+xKbbL6zcXM73pYicGPGUKKVUuMXDQ3qUijNuWvwTgB+IyA6srn4BjDFmdERTppRS4aItxsSm2y+s3AT+iyKeCqWUioRYB4xYf39HofkYVm6ezrcL62Y659j/V7lZTiml4kasAocGLBWH3Nyr/2Hgl8CD9qR04G+RTJRSSoVFrK/j18AfHpqPYeWm5X4lcDlQCWCM2QN0i2SilEoKPh8sWRLrVCQHDfxKNXAT+OuMMQYwACKSFdkkKZUkXn8d7r4bvorpM686tlgH3lh/f0eh+RhWbgL/OyLyMtBTRH4IzAP+X2STpVQSyMuzXouLY5sOFTkasFQcanVUvzHmKRE5DygDjgemGWM+jXjKOipj9NpiZdFyED3a1Z/YNB/Dys3lfNiBXoN9OGjgVyp6dHBfx6D5GFYtdvWLyCL7tVxEyhx/5SJSFr0kdjBagJWKnljvb7H+fhUZ8+fDnXfGOhVtFuyxvBPsVx3BH05aESi/jRtjnYLkoS3+xBZv+Xj//dZrgvbgthj4RaR3sAWNMSXhT04SiLcCrGJj927Yvj3WqUgeGvhVJHS0wI/1GF6DdW/+/sBB+/+eQB4wKNKJ65C0IlAAhw7FOgXJIdaVsu7vHVuCbt8Wz/EbYwYZYwZjXb53mTHmCGNMNnApMDdaCVSqQ4p1QEoWsa6Yfb7Yfn9HEevt2JIE3b5uruM/zRgz2//GGPMx8O3IJamDS9CColRCi9fAodyJ1+2XoPW5m8C/R0R+LSID7b+HgD2RTliHFa8FWAVXWgrTp4dvR9dyEF16jl9FQgcO/JOBPsC/7b8j7WmqLbQiSExPPmndYnfRolinRIVCr+PvGOI1H+M1Xa1wc+e+EuCnUUiLUvGrttZ69Xpjmw4VmlhXzAnaIow7sd6OLUnQ7dtq4BeRPsD9wEgg0z/dGHNOBNPVccVrAVbRpeUgujS/VSQkaOB309X/d2AT1uV7jwA7gWURTFPHphWQUslD9/eOLUG3r5vAn22MeQ2oN8b81xhzK6Ct/bZK0CNEFWYd4XK+L76A6upYp8KdWFXQur+HR7wFWP/+m6Db103gr7df94rIJSIyFgh6Vz+lOpxwVzxtWV9FRfwE2h074Gc/g9//PtYpCU4H96lISILA/5iI9ADuA34BvAr8LKKp6si0IoisPXsgJwc+7YAPkzzrLLj44linwlJZab3u3h3bdLRG97eOIV63Y0cM/CKSCgwzxpQaY9YZY842xow3xsyKUvo6nngtwB3F5s3W65w54V1vvHTNl5fHOgWJSVv8iS1e8zFe09WKoIHfGONFr9kPrwQtKEolND3Hn9jitd5M0O3b6uV8wGIR+TMwE6j0TzTGrIxYqjqyeC3ASnVEse6p0f29Y0rwc/xuAv8Y+/VRxzSDjuxvG60IIkvzN3mUlcHtt1t3VRw4MPA8Wh5UJCVo+XJz576zo5GQpJGgBSXhxLqlpyJv0SLYvh3+8hd45JHg82pXf2KL13ozQbevmzv3/TzA5FJghTFmdSvLXgg8B6QCrxpjHm/yeQbwV2A8UAxca4zZaX/2IHAb4AXuMcbMsae/jvVo4APGmFGtpT/uxGsB7mji4fK7aK5PBaeD+1Q4JXhXv5vL+XKAu4Bj7b87gQuB/yci97e0kH1FwAvARcAIYLKIjGgy223AQWPMUOBZ4Al72RHAdVi3Cb4QmG6vD+ANe1pi0opAQcJWGHHFzb6k1/F3DPGWj0kQ+PsB44wx9xlj7sNqnR8JTAR+EGS5U4BcY8x2Y0wd8DZwRZN5rgBm2P//CzhXRMSe/rYxptYYswPItdeHMWYhUOLmx6kkFKku/nCvN94qso5K87ljiNftGK/paoWbwH8kUOt4Xw8cZYypbjK9qWMB59098u1pAecxxniwTiFku1w2MSXoEaIKMy0H0aXn+FU4JXiL382o/r8DS0XkA/v9ZcA/RCQL2BCxlLWTiNwB3AHQv3//GKfGIUGPEFWYJWiFoUKk+3t4xGs+Juh+3GqL3xjzO6wAesj+u8sY86gxptIYc32QRQuA4xzv+9nTAs4jImlAD6xBfm6WbS3drxhjcowxOX369All0chqrQCXlsJvfxs/92RXqiPQc/wqEhJ0+7rp6scYs9wY85z9t9zlupcBw0RkkIh0whqs1/RWv7OAm+3/JwHzjTHGnn6diGSIyCBgGPC1y+9NbC+9BB9+CP/5T6xToiIpQVsKKkQJGhiUSwm6H7sK/G1hn7O/G5gDbATeMcasF5FHReRye7bXgGwRyQV+DjxgL7seeAfrVMInwBT79sGIyFvAV8DxIpIvIrdF6jdERGsVQUevKObPtx6ic+BArFMSmnBvl0hVGMXF8ItfND5EJ1ISrZxqiz+xxWs+Jmjgd3OOv82MMbOB2U2mTXP8XwNc3cKyvweaPfPTGJPYzw5oraB09BvPvP++9bp1Kxx5ZEyTElPOiiycldqrr8KCBfDRR3DNNeFbb1PxVBG72Wc08Ce2eMvHBB/cF7EWv2qB2xZ/Rz8ASHbOCiPeKjU34qHC8+ebm/zTwK/CKdb3h2gnDfzxSgN/fAn3jh6pFn+0xEPgTwSJuG3jUbzmY4LuBxr4oy1eC3C0RSofIp2/1dUwd27715PogT8e0uw/GIvng+R4yKeOIN7yMcG7+iN6jl8FEG8FONqiVUlH6nueegrKy+Goo+Ckk9q+HmeFkYiVRzykWbv6Vawl6PbVFn+0Jfuo/kRXXm69tnfUfKL2eET7e8JFA7+KBK831iloEw380ea2Iojn7stkFMnL+eKh9RyqeEpzPO8rGvjDI97yUQf3qZAkaEFJGImSvzqqP7q0xW95910oCOkmqPEh3vLRL9H2A5ue44+2ZO/qj/SRcqLsiJEa3Bet1m+ilVMN/FBXB3/4g3X/jNmzW59ftS6etm8ItMUfbfHU1e/zRT9QRvp3RWpHDHe6tcXffomQb/GURv/56EOHYpqMNomnfISEH9WvgT/a4qkAT5oEEyfG5rsTvcXf3gMB5+93k+Z4KjcQf+lpjbb4I7NvvPACvPFG+NfbVDzlo5MGfuVKPBXgvDyoqYl1KsIrUXbEULv63Zab9pavujrYuLH1+RIln2Mtnvb3SIxA/8tf4M9/Dv96E0U8bd8QaOCPtgQtKAkjlGu7YynUrv5QA21beyR+/3u48UYoLAxveiIhlG0cq/IQD/nkl6CXnsUl//6VoHmqgT/a4j0gRUukKsRIrbfpdot2V3+ov6ut5WztWuu1tfsUxENACyUNut8lbJAC4m/76eV8KiQJWlDCJtKDYhLlBjahtvjjaVAoNKY/ltfQJ0Lgd/O9W7fCvHmRT0siB/54FQ8HwG2ggT/eJMuBQaQqoWhVbu3d4QMF/spKqKiIzPcFU18Pf/xj410J3YiHp0gmQqXrZn+ePBkeeCDyaUnkwB9v9WKCt/j1Ov5oc1tZJWiBci3SLf5IB6RwVqL+NH/nO9br8uXN54lkufn4Y/jb36yBfW4rtHgIuvGQhtZEKo233gpnngm33OJ+mQQM/GXl5VRWVnJMvNaHbd2+BQXQpw906hTe9LikLf5oc1uA43EnffBBq2XYkkOHYMkSd+uK1O+LVEXb9EAinC3+SJzjD+XAx+OxXuvrG6clwo2mQknDo49Gp1UdLWvWWJfShSIe65RWbNmyhYJ4vtNgW+qBqiq44gprIG2MaOCPNreVVTy2Zj791GoZtuSee+Duu62WY2va+/s++ADWrQv/et0KZ+B3I9RA25bA7FzGfzDQkkDpj3aZ9Qcytwc50TiP3lQ8HCD5xWOd4lY85eMNN0BZmfV/W9JVWmq9LlsWvjSFSAN/vPEXpETcSXNzrddgQSNcg/t+9zv4wQ+aT49WBdFaYGxNvI7q92+fUAP/3r1wyinw0Udt+962iPdLN0tK4OuvY52KRu0ts8qyaVPj/23pRfGPpcnKCk962kADf7QlcovfLTct/kTr6m+63SIxuM/t/MGE6/ri1oJE0zTv2GG9zpnTvu8NRbzvI7fdBjNmuJ8/0r8nAbv6415bDjr9vQVduoQ3LSHQwB9toQT+nByYNq393+nxWE/lilZF6TxX3JJIVULRav21Ny8jfee+UNIXaEBfqC3+YK3vjz6CBQvcp8ctfxpKS+H734fdu8P/He0RanrctsjbWvbiIfCXl1vjE0Lki4e0B9KWdGngT0Kt7bRNK+5wPEXrH/+wnsr1/vvtX1cw/gBy0UXw5JPB521P4AwWBKN1cNPeiqi1Fn9dXWMF0XT+UNffGuc5crc9BqEcYD38MPziF+7nd8v/G7/4ArZsgTffDP93RJPbwN/WLvtY9pDk5FhjgH76U+uKhBB/g8dNYyIW3PRuNuV/SFLnzmFNSig08MercO6k/sEk/tdoeOedwNPD0RUdbNlECfytneO/5x4455zA87vRlnxoT4vfv12jeV1/S70OiSSUPA91vnAtFy5ffgnr11v/hxgw6ztS4K+qsl5jdCkfaOCPPreVUzh30tauzY5FS8D5nXv3wqpV7pcNtrOFq/K/6SY477zG9+G+nK+1rv6m1/JHq8Xv15ZR/dEWzTS89x786ldtWnTf/v2Ut3RjJmcvhdt9vq1BMJzd5ZWV7etBDPE3eGN90NKStgR+/zIpsQu/Gvijze310dF8al40HwsbqMV/xRXwwx+6X0ewSiBcwWDDBjh4sOXPI93V3/SzUAcDhpI+5/ra2tUfar6//TacfTb87/8ePr24GO6/v+U7GAZLgzGwcGHwtLT2u+rqAi//v/8Lc+e2nqYA8vPz2bx5c+AP589v/D/SLX7/b3dTB7W2Pf/wB3jssdDT4P/uEANmVLr6q6qgqCi0ZdoS+Gtrrde5c+Gtt0JfPgw08Eeb2zuiVVeH7ztb6351EyTcVOyBvqeuLvBvDvUGNk3X2ZJIdfdGclR/sC7r+nqrLCxa1HJaAi0bSvqcN/DZvv3waS1pun7//G7yv7wcnnrKen3vPZg+vbEyfO01Kxi6uSywabn97DP4+c+DV6at7Vff/nZ4BtQ6L/kKxrnPtLXFn59vPWK7NW7LxB//aF2a6Z9/xgyYOPHweeynNxYUFLBv377g6wtUJsrK4MorXff0RaWr/4Yb4MILQ1vGX27d2LYN9uw5vP56+unQvi9MNPBHWywDf0vdy24Cf1uObAsLrYr0X/9q/ll7bgATLC3RupwvnOf4gx1U1NRYrSvnYMlgv7E9gf/jjxunhdrib2mbNJ0vL6/5MwFef71xTEgoY0Cartu/Xsed3srKyti5axcNc/rPr779NqxcGXh9n3zi/jsDKSiwggjgC+VAtK0t/u9+F666yv1yrTUE/v5369Uf1J5/3sq3AGVq77595Ld2Zz3ntvR/94YN1lUPzz3XeroBT1t7ObZvh+99z934Jv/BU1mZdVrHPwgvmFAOSK69Fi6/PLSDhQjRwB9twSqCLVvAf/TsDPx794b/u0MdVBTK7Vz9/BWCM6AEq9hbO7gwxqqsg+1sPp9V2X7+edvPQfpvRBRMGLr6PV4v+QUFzc9fOn9fXR3s3Nls2VbT1ZbA39q0YGloKag411NfbwWo++5rvj7/qa3U1MDrd5MGv/37G/5dtWoVRUVFjYHDH/ifegruuOPw5dxUyG72FccpIhPKdnBbptwGmw8+CO1grqmmpxuDnX4MVicE2q/901zWJW3u6n/tNdi1C776yv0yf/2r1Q3/73+3nB6vl+UrVrC5DZcmauDvSK66yiowLdi5axerVq2Cls71gXUtsr8V4gz8f/lL8O/Oy7Mul2nagtmzBy6+uPHAoaVg76ZCaFqBt5W/IgxUIbbWy/Huu1ZlHeRc64F9+1i5ciXFJSVtOwcJ1o1XnL75xjp37BSGwX3+btJvVq8+fH3O/A1USTi345Yth8/jX08olXw4An9LZcKZtspK63Xr1ubz+X+Tf8BTewL/f//b0FozDbPa87bQcgXc9bK5qbQdgc7X2u9wHii9+y48/vjhl3EG4nb/+93v4De/aXzv9hy/X9P8CJY/oQb+EINfa4H/rQcfpPDdd0NaZ4v27LFee/VqcZZq+wDyq//+N/T1t6X3NMw08IeDz2cF3z/9qcVZioqK8Pp81jyOFkmD/PzD3zuPrrOzg3+/fwS48+ge4D//gQMHqP/Pf9i+fTuVzgFTXi9eu9XpKkg4dzy3FU+gAh4sMDVtURQUWF2Cfv6uOP9d4gIosG+acjDYwLzWOCs4n6/5gQC0v8Xv9TYEBU9dXcsHVoEqSJ/PCg7bt1sHi85z0k1b/IECT0FBY8sXAgf5Vn5ffV0dubm5lPjzuaXKLFDgDyScgR+aBc+GG8BUVbX8+OGmgc3rbRzz4Odm0K3ju0Pq6n/7beu02EsvBZ8v2EHZwYOH7zNOoR6sNv2tzjITSpocZSM/P5/lK1bgC3H/DNbVv3v3boY9/jhFN93U+opqaw8fUBmIi4cC+bdrZgiXr1ZVV1NTW6st/g6jlXNBpY7zSz6fr/n5ptJS6zwd1mMol69YwSFn936wgjJnjtXCgebdrPb7/fv3U3LwIIs+/7zxM4+HNWvWsHr16uCV/Pbt1g7vDEb+71u50nooT0vL+ytSR+VnPB42btrE14G63vzzz5ljPXnsiiusy+r8/EEhlC5HN6PDg2kp790OiGypsnUEB2PM4fnrrORaavFfdx1cc4313nlXPGfALyqyBmnNmnX4sldcYd0nwC/QgVwrLf5dO3ZwqLSUxYsXH76OpvnvNvD70+3fxm661IMFsibb3ets8bd0vrdp4P/rX6083rKlcZr/9+zdCxs3Bl6PM/Db31sH8MtfurulcaD0lZdb371xIzzzTODl/AepN90UuAVu56nH62W2mxuDvfACLF3a+L5J/hz2DcH2B0cZ2GPXa9UHDrT+/Q7BWvw1doOg3uOx7hOQk9P8ron+7f/cc9ZVI8G66P0HJUEOdPwHIr39d1hdvNjK8717rVMLU6ceXk95vWzYsIF169Y1P0g+5RR3gzPDSAN/OBQXN/4f4Dad9/zkJw3/ezye5hVgSYnjX+v/um++afzcX5EUFjYPBA89ZBU6pxYqRHEWOK8Xr39nbWmn9Xqtiu/nP28Ye2AAM22aVch/9SvrMbyBxiB06oSpqODAgQPU1NRYvRyPPkplURGVlZW8G+gpf/4d5aGHAp/e8AcFf1ccWN3+W7daD0N58EEGOA9uAM46y31Lp6zM2omd8595ZuB53azz3nvh/PMDflRTUECxXW5aa/EbEXzGWHnv/25nxen1Wq1F52VYPl9ja/U//2mc11/2Vq9unNZSV/+WLY3nRv3rqquDPXsaz137A0xLgdpZ+dnB2ABbt26lzNnyzs21Ro7v2mW9D9a6tPm8Xg4UFjYGdacmgb+hy726uuWu9Kbjavw9S8688u9/l10GN94YeD0BAn8KWFcdPPRQ4GVaSsfGjdalhGefbQX1H/3o8PQ49+nKysYA4twn6+pg7dqG/XxXbi5/u+QStm/bdvj3GgP//Gfj+//+F6ZMaZ6u11+HFSsOH78QLPAHCNq1gXo9S0paPGXgq6xsLNvvvdfY+PjiC465/fbG9Pu7+/3l1t5vWLHCeu8fL+PfRuXlVk+p83v9B16ByqC/l84u7wP9+9OMGdY+eNllcNdd1hijCRMaTxE6t6m9zbz+8Ug+n3X3ySiKaOAXkQtFZLOI5IpIs4dhi0iGiMy0P18qIgMdnz1oT98sIhe4XWck+bxePr/4Yr66/37qa2p4b8IEdixeDCUljUe/V1552DKmpoZtjpt0eLxeq9AtX26d7//rXw+/XtwugIdKSxsrtPJy8HjYP24ce2+9tXFer/fwrkQRKCrCd/LJmA8+aDh94O8HSLF3wPpp0/D4dxbAOHbMO264ga+nTrUCoH/MwvLlcNdd+Ixhx44drFy50joI8Qfit9+2v8B67/X5MOnpbF2zhrzdu1myZIl1imPWLOrtp5UNse9aVVhY2NiN11JL3l/Z+ns0cnM5dOgQ1TU11sHH5Mnw4x/Dp582DJSrqqqy8howTz7JY9/7HotnzwaPhwNz51oVwYYNVn4bY+2YQe6tXl1Tw85du6iqqiJ32zZq/ZVCTc1hB14Vq1ax8cYbrYrwyy+t3iCP57CKZfn8+ex0PCK2qrISPB68Pp+1zZ2XONXWsnXLFlauXMk333zD1q1bAx90PPUUFBVRW11NwZ49eP7+94ZLrvB6G3+b4yCz4QDK42koR7W1tezZswdfXR1VV15J6c03WwHrmmusSuzpp+Hyy6m1DzyG+C819JchR6+TMYYKfyCHhp6xmpoaSsvKyHUMovR+/jn7tm/H5++GdRwcV1ZWBuzq3bBuHXl5eexxHgj6NenObzjIraw8rIfuy0WL+I3/dsJ25bx//372n3oqdO1qTV+7tvGgy1FGfcaw1d/q37sXz223Wb0s/nynMfCnORPzwAPBz5c7ryWfNs0KdGBdIlhRQXVNDZX+/Pn1r6mpqWHvvn3UO3se7R5EAJ59Fm65BXbtoqKyksrKSn4O5DV9TPHOnfDEE1RUVLArL4/9+/cffnA2b551w6Hp0wEOO+Dav2MHbz3wAGVlZVTZ+8aypUuprKggf8WKZqc8qvftY2tubkNDh7w86yB5+nTrgNPZywL0e/FFq4yDdSB0333W/rtiRcM+b6Bx2z75JPzgB2zNzWXFihXUzJwJ1dXWQTZYwf3ZZ60Dqt/8BjZupLaujorKysayt2GDdVWS/6C3srJh29bb62kIoCtXNl6a59j+fPQRZaWlfPTiiw2TahYuZPfu3axatYqNGzZQW1cX9bv4iYnQdc8ikgpsAc4D8oFlwGRjzAbHPD8GRhtj7hKR64ArjTHXisgI4C3gFKAvMA/4lr1Y0HUGkpOTY5Y3vRNaG83u0YPeFRX0GjiQcrtAlF9yCb3nzOG4/v3p1KkTXe+6C7p3p3TiRLbk5CCOHXLAgAGkTpzI7rfeom/fvmR27kzhqacycN06aycpKGjcqQFz9NH0FaFnz55stCuZ8du3c2jhQnI/+4yUN9+kU0YGAwcMQFJS2LlrFzXV1Rx77LGICCJCbU0NBwoLKTn1VMa98AJFZ51FhaNFNODll8kYPBhfSQmrrr2WbsC4ceMoKS4mIzOT9LQ09u7daw2Ys/UZMYK0gwfp3asXdfX1dHv/fbw33URpURF5eXl07d2bg9deS+r06RQfcQQX3HwzZsEC9hQUsHffPsq6dCH3xhspf/llzklNZeyYMRzo2pUdvXszdPVqunfrRmlZGd27daPTnDl885OfMGTdOjI7d8bn87FmzRpSRBh+wglU25XowYMHSU9Pp9Cx4w0ZMgQRITc3l+0pKYweP56KZcs4IjubrKwsKsePp8/GjXTJysJz5514nnsOr9dLp06d8NTXU1pWRo/u3dnQpFu359VX0+fMM+m2ciWydi37amspvuwyMqZPp7S0lJ4vv8yQV17BGEN1dTVrDhyg9KyzOLqkBE+TXondp5/OqPvuo3TSJATI6toVYwzpaWl06t6dA00OSGqnTmX0rFns2rWLoUOGkJmZCYBvyhS+/MMfyLQHz/XOyeGY6mrq6+spLCyk/v77GX7BBZRecw3FJSV4PR7M97/PsQsWsGnzZpx1wdH33MM+e8zKkTfcQPaKFXgvvphuW7ciBQUcOHCAPEe6evXqRd++fek8cSI7e/fGt2IF73k8nDVvHtnZ2XTr1g1vair78vLo0qVLw+mvIUOG0LNnTzZt2kRlZSV9+/al7zHHWCu96irqr7+eNwcOpOvYsVy9ZAklX35J5gkncGD2bLIefrghDcOGDSOrSxfq6uqorqmh9xFHUG0MG+z9Pi0tjfS0NPpNm8a2xYvpvmgR/Y47juc3buRsr5ey++7j7PPPR371K5bbrcP+l1xC1ZIl9O3bl/yCAkpKSkg94QQ6DxvGsZs3U1JczJdFRYydMQOZMoXyigr6jRvHEUOGULN2LVVVVeQ5unDHjx+PYB0MFB1zDOmlpfRKTaUkPx+MoXfv3hSXlJC/Zw9DZ81i5Usv0X/RItLT06mrq6Nv375UV1dbB3/A4MGD6dmzJ2vXrqW+vp6USy+l07x5dEpPZ8jQoaTaB+KmTx/qCgrw9O/PRsf9EfIuuYQx//d/9M3O5k/XXsvlI0cy5Isv+KZJN3j//v0pKy2lpqaGmtpahgweTM9evSgtLW04eFuSlcVplZV8CEhGBjdcdBEH338fb48epNrbuu8xxzR09R911FHs37+frC5dGPzpp+x9/32qXniBo446qmF7FfTtS5Vj8NyIESNYuWsXfTMy6JqVRVbXruzfv5/09HT22z0IgwcPJiUlhR49eiDAtu3bOXjwIH379iXtppvIe/xxevTowRF9+tAtK4ua2lp27NiBiFg9k0DXrl2pqKjghBNOIKtLF+ugacoU1g8YgG//fvr06UO+3ajq0rkz/QcMICsrCwEqKitJT09HRCgvL+fQoUN4PB7KWxpXYjv+Jz8h48kn6WTvy+EgIiuMMTkBP4tg4D8d+K0x5gL7/YMAxpg/OOaZY8/zlYikAfuAPsADznn989mLBV1nIOEM/HPPPJPezpupBJBiB1znEfHuY47huDZclvcVcHrISwVWBnRv5zr+B+sorCWzsY7QhjaZ3jUry2oJiAQc6dytW7eAO0dqairpaWnWoJgY+g9wWYDpKSkprY/cjrDu3brh9XqpdNE9npaW1vZrol1IS01t6Glpq8zMTFJSUqitrW1sqbu0EhgH1kFdsDEFAaSnp0f8RjGdOnWiztE93yk9nTr7O1NTU1393qeBC4GRAT4rAXrb/6empFiniCJQx2d06mS1VONUWloaneztWW+X97bsq5kZGaSmpZHfrRu9mtyoaCswrI3pywUGcnhPUMnkyZz/j3+0cY3NBQv8aYEmhsmxgLOpkg+c2tI8xhiPiJQC2fb0JU2WPdb+v7V1AiAidwB3gHXEGi6nPfEEi378Y/auX08vj4dsYE1GBhX9+/Od3bvZOnQoPXr2JGvLFrqWlTGvd2/u/OwzcgYPZstll1E2dy47Ro5E9u5ldkkJZGXxw5NPpsvixcw58kj2HX00x48cSfbMmTxXW0shVpfHFuCgnTnnAnuBb3XvTuaddzJt2zaOrK+nq8fDoK1bGWUMhVlZHLz+etbm5XHOrFkcUVhIqgiZxnAAqBoxgk/37+ecPXs41lExfIB19HVPRga1tbVsAYanpPCJz8fGbt0oLS/nfuDifv1Yt38/GfX1fAvrdEJPIO/yyxkwcSKLn3mGsZWV7Dj+eNbv3cu5RUW8168fOfv30+O88/DNn0+Xujr6er2kde/OF6mpDBs0iPySEpaWlnJdZiadamowxhzWTViTnc38sjLSsrM5WFLC7kGD2LFrF/VeL4/deCPZs2ezeN8+CgYO5Mg9ezi3ro4e3bvzROfOdN+/nwXA5CeeYM+vf81F9fUU3nIL2atXc4Tdvb7ymGM4PiWFLd26UVhRwdayMq7p3JnBl1zCRq+X6jlz+HdGBjcXFzOitpZ9Ph+9gYKuXdmXnc3So4/m56NGseeTT5hRUMBVwABgtwhb0tI44PEw3hheSUnh2/ffz49ycth3002k1dbSuVcvfjV0KD3WreOa669n/fbt9O7dm7MuvZTCf/+bPWvXsnXfPsanpFA5dCgfrliBJyODzSL8tKqKwcB84KKMDDI6dWJjdTXf9OjBiSecwLG9enHkV1/h9XhYcfTRLDrjDAo/+ogLKiupKy+n/+DBHJOXx6Zvf5vtPh85eXnUjRqF7+OPqTeGT4HzBgwga9cuCoAjsLo5a4GNWAeDx2Vmcm1GBl+WlnLkhAmY5ctZ1aMH37v2Wo5/800KjaFn587sGzCAsT/6ETuzsqi54Qa+qatjps/HzUAvoDY1ldQuXTjGGHoZw8tDhtAzI4NRK1cyBNickoLH5+NbwILOnbm8Vy98xcX4vF4y77yT1bNnMz4vj/quXdlz5JEM2L6d1NRUVh9zDMfv2cPirCy+HjuWSevXYw4eZEyXLtTW1FCRns6K0aPpu2oVg42h2BiWAielpVHp81kDuYYPZ2F+Pt+tqKC0Sxd6VFVhMjLo+a1vUXv00Xjmz8f06cPOqVOpzsvj+Fdf5RFgSF0dV9TXsy81lTy7TCxJTeXaTp2Y3rs3lQcOcGaXLmRkZkJGBql79tDX5yMTmNOrFyY9nc7XXMOgZctYtGMHQ26+ma5PP80An493hgxh0IQJ7Jsxg9nDh3NnXR1jtm/H6/OxNyWFzKws9qekcKC6moxJk+ickkJtSQmD162jtKqKrjU19E1NpbMxlIhQ5/FQVFXFauCkyy6jat06fN/7Hj1zcii57jpGA8tSU5nQrRtrevakJiUF38GD7CkrI7tLF8alppJeW8tNdXXc0asXByorOTY1lRFAUXU1xwLSuzfp5eV07tyZuuxsyioqqL/gAp6cOZO+Z5zB1Ucdhad7d0pmz8ablkZVSQmdu3ShtLaWzYcOUQz8unNnPjaGE2pq+BooOPpoxtfUsKZvX87av5+xtbUNB1X/7dyZiamprM/OZtHppzN+3jx6VlUx0uOhYvhwyn75S3Y8/zxZS5dSkpND3caNnFNRwfLaWo7u0QNPdfVhjbA5WVn0njKFPz39NBd6vQ3d0WA1YrbU1LBh9GhSKiro0aMHA9evZ+vw4WRPnEhtv36ceOaZpO/YwexbbqE6I4OTBw9mZCSeYNmCSLb4JwEXGmNut9/fCJxqjLnbMc86e558+/02rED+W2CJMeZv9vTXAP+1akHXGUg4W/zh5PF4SEuzjr2MMUgrl4bU19eTnp7eOMGYsD4NzefzkdLKgyMCpdPr9TacVmjtN7SJ/3faR+vVtbV07tz5sPxrljdt4PV6SfXfQKbVJLW+vVpSV1dHJ8c5PW9dHeLzkdKGbj5/OozHY43lSAtyLB+gvPj3/4hst6a8XmsciPO77DT5y159fT0pKSmutkPD9jIGr8dDSmoqEqT8er1ejDENZaY9asrKyOze3XX+OctLbW0t9fX1dO3a1SrTKSkhlT2/6qoqUkTICPB4V+e+0RbByrfxepEQ0+pXV1fX0BUeTu3ZH1vi8XhISUkJWifW19eTlpbW8N01NTUNp95iLVYt/gLgOMf7fva0QPPk2139PYDiVpZtbZ0Jw7ljuim0zQJbmAt6a0Hf+srm3xlqhRUy/3fa6etsV3TO/Gtv0IfQfkd7KplOTQbypLZjYI8/HeKmkg+Q5qgEfL9A+Wt/v7/shbIdG7aXCKkulgtnOc3s3t3+anf555wvIyODjIwM6439u9uSts5durT4WXsPboL9rrYGfWhe9sMlEuXYTR42La/xEvRbE8lR/cuAYSIySEQ6AdcBs5rMMwu42f5/EjDfWIfQs4Dr7FH/g7BOpXztcp1KKaWUakHEWvz2Ofu7gTlAKvC6MWa9iDwKLDfGzAJeA94UkVyscSnX2cuuF5F3gA2AB5hijPECBFpnpH6DUkop1dFE7Bx/PInXc/xKKaVUJAQ7x6937lNKKaWSiAZ+pZRSKokkRVe/iBQCu1qd0b0jgKJW51KBaN61j+Zf22netZ3mXfvEIv8GGGP6BPogKQJ/uInI8pbOnajgNO/aR/Ov7TTv2k7zrn3iLf+0q18ppZRKIhr4lVJKqSSigb9tXol1AhKY5l37aP61neZd22netU9c5Z+e41dKKaWSiLb4lVJKqSSigT8EInKhiGwWkVwReSDW6YlHInKciHwuIhtEZL2I/NSe3ltEPhWRrfZrL3u6iMif7DxdIyLjYvsLYk9EUkVklYh8aL8fJCJL7TyaaT+nAvtZFjPt6UtFZGBMEx5jItJTRP4lIptEZKOInK7lzh0R+Zm9v64TkbdEJFPLXctE5HUROWA/YdY/LeSyJiI32/NvFZGbA31XJGjgd0lEUoEXgIuAEcBkERkR21TFJQ9wnzFmBHAaMMXOpweAz4wxw4DP7Pdg5ecw++8O4MXoJznu/BTrMfd+TwDPGmOGAgeB2+zptwEH7enP2vMls+eAT4wxw4GTsPJQy10rRORY4B4gxxgzCus5KNeh5S6YN4ALm0wLqayJSG/gYaxH0Z8CPOw/WIg0DfzunQLkGmO2G2PqgLeBK2KcprhjjNlrjFlp/1+OVfkei5VXM+zZZgDftf+/AvirsSwBeorIMdFNdfwQkX7AJcCr9nsBzgH+Zc/SNO/8efov4FyJ6nN244eI9AAmYj34C2NMnTHmEFru3EoDOov1ePQuwF603LXIGLMQ68FyTqGWtQuAT40xJcaYg8CnND+YiAgN/O4dC+x2vM+3p6kW2F2AY4GlwFHGmL32R/uAo+z/NV8P90fgfsBnv88GDhljPPZ7Z/405J39eak9fzIaBBQCf7FPk7wqIllouWuVMaYAeArIwwr4pcAKtNyFKtSyFrMyqIFfRYSIdAXeBe41xpQ5PzPWpSR6OUkTInIpcMAYsyLWaUlAacA44EVjzFigksauVkDLXUvs7uUrsA6e+gJZRKnl2VHFe1nTwO9eAXCc430/e5pqQkTSsYL+340x79mT9/u7Uu3XA/Z0zddGZwCXi8hOrFNJ52Cdt+5pd8HC4fnTkHf25z2A4mgmOI7kA/nGmKX2+39hHQhouWvd/wA7jDGFxph64D2ssqjlLjShlrWYlUEN/O4tA4bZI107YQ1+mRXjNMUd+1zfa8BGY8wzjo9mAf5RqzcDHzim32SPfD0NKHV0lyUVY8yDxph+xpiBWOVrvjHmeuBzYJI9W9O88+fpJHv+uG1lRJIxZh+wW0SOtyedC2xAy50becBpItLF3n/9eaflLjShlrU5wPki0svudTnfnhZ5xhj9c/kHXAxsAbYBD8U6PfH4B0zA6uJaA6y2/y7GOgf4GbAVmAf0tucXrKsltgFrsUYWx/x3xPoPOAv40P5/MPA1kAv8E8iwp2fa73PtzwfHOt0xzrMxwHK77L0P9NJy5zrvHgE2AeuAN4EMLXdB8+strPEQ9Vi9Tbe1pawBt9r5mAvcEq306537lFJKqSSiXf1KKaVUEtHAr5RSSiURDfxKKaVUEtHAr5RSSiURDfxKKaVUEtHAr5SKCRFZICI5sU6HUslGA79SSimVRDTwK6UaiEiWiHwkIt/Yz2a/VkSmicgy+/0r/iex2S32Z0VkuYhsFJGTReQ9+9nij9nzDBSRTSLyd3uef4lIlwDfe76IfCUiK0Xkn/azHhCRx0Vkg/0c86eimxtKdUwa+JVSThcCe4wxJxnr2eyfAH82xpxsv+8MXOqYv84YkwO8hHWL0inAKOAHIuJ/YtvxwHRjzAlAGfBj5xeKyBHAr4H/McaMw7r73s/t5a8ERhpjRgOPReYnK5VcNPArpZzWAueJyBMicqYxphQ4W0SWisharAcHjXTMP8ux3HpjzF5jTC2wncYHkOw2xiy2//8b1m2dnU4DRgCLRWQ11n3OB2A97rUGeE1ErgKqwvlDlUpWaa3PopRKFsaYLSIyDuv5Co+JyGdYrfgcY8xuEfkt1r3a/WrtV5/jf/97f/3S9L7gTd8L8KkxZnLT9IjIKVgPjZkE3I114KGUagdt8SulGohIX6DKGPM34P+wHm0LUGSfd5/U4sIt6y8ip9v/fx9Y1OTzJcAZIjLUTkOWiHzL/r4expjZwM+Ak9rw3UqpJrTFr5RyOhH4PxHxYT157EfAd7Ge2rYP6/HUodoMTBGR17Ee9/qi80NjTKGI/AB4S0Qy7Mm/BsqBD0QkE6tX4Odt+G6lVBP6dD6lVMSIyECsxwuPinValFIW7epXSimlkoi2+JVSSqkkoi1+pZRSKolo4FdKKaWSiAZ+pZRSKolo4FdKKaWSiAZ+pZRSKolo4FdKKaWSyP8H2z93pjRdVrUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize=(8,3))\n",
    "ax = plt.subplot(111)\n",
    "ax.plot((grad1-grad_true).abs().sum(1), c='black')\n",
    "ax.plot((grad2-grad_true).abs().sum(1), c='red', alpha=0.8)\n",
    "ax.set_ylabel('gradient error')\n",
    "ax.set_xlabel('samples')\n",
    "plt.legend(['torchdyn', 'torchdiffeq'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Integral loss (backsolve)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 302,
   "metadata": {},
   "outputs": [],
   "source": [
    "t_span = torch.linspace(0, 2, 500)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 303,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "tensor(8.4026e-06)"
      ]
     },
     "execution_count": 303,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from torchdyn.core import NeuralODE\n",
    "\n",
    "def reg_term(t, x):\n",
    "    return 0.1*x[:,:1]\n",
    "\n",
    "x = torch.randn(1024, 2, requires_grad=True)\n",
    "\n",
    "# solve with autograd augmentation\n",
    "node = NeuralODE(f, sensitivity='autograd', solver='rk4', interpolator='4th', \n",
    "                  atol=1e-4, rtol=1e-4, integral_loss=reg_term)\n",
    "x0 = torch.cat([torch.zeros(x.shape[0], 1), x], 1)\n",
    "\n",
    "t_eval, sol_torchdyn = node(x0, t_span)\n",
    "loss = sol_torchdyn[-1,:,0].sum()\n",
    "loss.backward()\n",
    "grad_autograd = x.grad\n",
    "x.grad = 0*x.grad\n",
    "\n",
    "# solve with integral loss without aug\n",
    "node = NeuralODE(f, sensitivity='adjoint', solver='dopri5', interpolator='4th', \n",
    "                  atol=1e-4, rtol=1e-4, atol_adjoint=1e-5, rtol_adjoint=1e-5, integral_loss=reg_term)\n",
    "\n",
    "t_eval, sol_torchdyn = node(x, t_span)\n",
    "(0.*sol_torchdyn[-1].sum()).backward()\n",
    "grad_adj = x.grad\n",
    "\n",
    "(grad_autograd-grad_adj).abs().mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 304,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "tensor(1.3780e-05)"
      ]
     },
     "execution_count": 304,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from torchdyn.core import NeuralODE\n",
    "\n",
    "def reg_term(t, x):\n",
    "    return x.sum(1, keepdim=True)\n",
    "\n",
    "x = torch.randn(1024, 2, requires_grad=True)\n",
    "\n",
    "# solve with autograd augmentation\n",
    "node = NeuralODE(f, sensitivity='autograd', solver='rk4', interpolator='4th', \n",
    "                  atol=1e-4, rtol=1e-4, integral_loss=reg_term)\n",
    "x0 = torch.cat([torch.zeros(x.shape[0], 1), x], 1)\n",
    "\n",
    "t_eval, sol_torchdyn = node(x0, t_span)\n",
    "loss = sol_torchdyn[-1,:,0].sum()\n",
    "loss.backward()\n",
    "grad_autograd = x.grad\n",
    "x.grad = 0*x.grad\n",
    "\n",
    "# solve with integral loss without aug\n",
    "node = NeuralODE(f, sensitivity='interpolated_adjoint', solver='dopri5', interpolator='4th', \n",
    "                  atol=1e-4, rtol=1e-4, atol_adjoint=1e-5, rtol_adjoint=1e-5, integral_loss=reg_term)\n",
    "\n",
    "t_eval, sol_torchdyn = node(x, t_span)\n",
    "(0.*sol_torchdyn[-1].sum()).backward()\n",
    "grad_adj = x.grad\n",
    "\n",
    "(grad_autograd-grad_adj).abs().mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "torchdyn",
   "language": "python",
   "name": "torchdyn"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
